diff --git a/lapack-netlib/SRC/zhb2st_kernels.f b/lapack-netlib/SRC/zhb2st_kernels.f
index a440b5c0d..2c0cb6870 100644
--- a/lapack-netlib/SRC/zhb2st_kernels.f
+++ b/lapack-netlib/SRC/zhb2st_kernels.f
@@ -1,26 +1,26 @@
*> \brief \b ZHB2ST_KERNELS
*
* @precisions fortran z -> s d c
-*
+*
* =========== DOCUMENTATION ===========
*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
-*> Download ZHB2ST_KERNELS + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
+*> Download ZHB2ST_KERNELS + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
*> [TXT]
-*> \endhtmlonly
+*> \endhtmlonly
*
* Definition:
* ===========
*
-* SUBROUTINE ZHB2ST_KERNELS( UPLO, WANTZ, TTYPE,
+* SUBROUTINE ZHB2ST_KERNELS( UPLO, WANTZ, TTYPE,
* ST, ED, SWEEP, N, NB, IB,
* A, LDA, V, TAU, LDVT, WORK)
*
@@ -32,9 +32,9 @@
* INTEGER TTYPE, ST, ED, SWEEP, N, NB, IB, LDA, LDVT
* ..
* .. Array Arguments ..
-* COMPLEX*16 A( LDA, * ), V( * ),
+* COMPLEX*16 A( LDA, * ), V( * ),
* TAU( * ), WORK( * )
-*
+*
*> \par Purpose:
* =============
*>
@@ -124,7 +124,7 @@
*> LDVT is INTEGER.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array. Workspace of size nb.
*> \endverbatim
@@ -147,7 +147,7 @@
*> http://doi.acm.org/10.1145/2063384.2063394
*>
*> A. Haidar, J. Kurzak, P. Luszczek, 2013.
-*> An improved parallel singular value algorithm and its implementation
+*> An improved parallel singular value algorithm and its implementation
*> for multicore hardware, In Proceedings of 2013 International Conference
*> for High Performance Computing, Networking, Storage and Analysis (SC '13).
*> Denver, Colorado, USA, 2013.
@@ -155,16 +155,16 @@
*> http://doi.acm.org/10.1145/2503210.2503292
*>
*> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
-*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
+*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
*> calculations based on fine-grained memory aware tasks.
*> International Journal of High Performance Computing Applications.
*> Volume 28 Issue 2, Pages 196-209, May 2014.
-*> http://hpc.sagepub.com/content/28/2/196
+*> http://hpc.sagepub.com/content/28/2/196
*>
*> \endverbatim
*>
* =====================================================================
- SUBROUTINE ZHB2ST_KERNELS( UPLO, WANTZ, TTYPE,
+ SUBROUTINE ZHB2ST_KERNELS( UPLO, WANTZ, TTYPE,
$ ST, ED, SWEEP, N, NB, IB,
$ A, LDA, V, TAU, LDVT, WORK)
*
@@ -181,7 +181,7 @@
INTEGER TTYPE, ST, ED, SWEEP, N, NB, IB, LDA, LDVT
* ..
* .. Array Arguments ..
- COMPLEX*16 A( LDA, * ), V( * ),
+ COMPLEX*16 A( LDA, * ), V( * ),
$ TAU( * ), WORK( * )
* ..
*
@@ -195,8 +195,8 @@
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, J1, J2, LM, LN, VPOS, TAUPOS,
- $ DPOS, OFDPOS, AJETER
- COMPLEX*16 CTMP
+ $ DPOS, OFDPOS, AJETER
+ COMPLEX*16 CTMP
* ..
* .. External Subroutines ..
EXTERNAL ZLARFG, ZLARFX, ZLARFY
@@ -209,7 +209,7 @@
* ..
* ..
* .. Executable Statements ..
-*
+*
AJETER = IB + LDVT
UPPER = LSAME( UPLO, 'U' )
@@ -240,10 +240,10 @@
V( VPOS ) = ONE
DO 10 I = 1, LM-1
V( VPOS+I ) = DCONJG( A( OFDPOS-I, ST+I ) )
- A( OFDPOS-I, ST+I ) = ZERO
+ A( OFDPOS-I, ST+I ) = ZERO
10 CONTINUE
CTMP = DCONJG( A( OFDPOS, ST ) )
- CALL ZLARFG( LM, CTMP, V( VPOS+1 ), 1,
+ CALL ZLARFG( LM, CTMP, V( VPOS+1 ), 1,
$ TAU( TAUPOS ) )
A( OFDPOS, ST ) = CTMP
*
@@ -281,14 +281,14 @@
*
V( VPOS ) = ONE
DO 30 I = 1, LM-1
- V( VPOS+I ) =
+ V( VPOS+I ) =
$ DCONJG( A( DPOS-NB-I, J1+I ) )
A( DPOS-NB-I, J1+I ) = ZERO
30 CONTINUE
CTMP = DCONJG( A( DPOS-NB, J1 ) )
CALL ZLARFG( LM, CTMP, V( VPOS+1 ), 1, TAU( TAUPOS ) )
A( DPOS-NB, J1 ) = CTMP
-*
+*
CALL ZLARFX( 'Right', LN-1, LM, V( VPOS ),
$ TAU( TAUPOS ),
$ A( DPOS-NB+1, J1 ), LDA-1, WORK)
@@ -296,9 +296,9 @@
ENDIF
*
* Lower case
-*
+*
ELSE
-*
+*
IF( WANTZ ) THEN
VPOS = MOD( SWEEP-1, 2 ) * N + ST
TAUPOS = MOD( SWEEP-1, 2 ) * N + ST
@@ -313,9 +313,9 @@
V( VPOS ) = ONE
DO 20 I = 1, LM-1
V( VPOS+I ) = A( OFDPOS+I, ST-1 )
- A( OFDPOS+I, ST-1 ) = ZERO
+ A( OFDPOS+I, ST-1 ) = ZERO
20 CONTINUE
- CALL ZLARFG( LM, A( OFDPOS, ST-1 ), V( VPOS+1 ), 1,
+ CALL ZLARFG( LM, A( OFDPOS, ST-1 ), V( VPOS+1 ), 1,
$ TAU( TAUPOS ) )
*
LM = ED - ST + 1
@@ -342,7 +342,7 @@
LM = J2-J1+1
*
IF( LM.GT.0) THEN
- CALL ZLARFX( 'Right', LM, LN, V( VPOS ),
+ CALL ZLARFX( 'Right', LM, LN, V( VPOS ),
$ TAU( TAUPOS ), A( DPOS+NB, ST ),
$ LDA-1, WORK)
*
@@ -359,13 +359,13 @@
V( VPOS+I ) = A( DPOS+NB+I, ST )
A( DPOS+NB+I, ST ) = ZERO
40 CONTINUE
- CALL ZLARFG( LM, A( DPOS+NB, ST ), V( VPOS+1 ), 1,
+ CALL ZLARFG( LM, A( DPOS+NB, ST ), V( VPOS+1 ), 1,
$ TAU( TAUPOS ) )
*
- CALL ZLARFX( 'Left', LM, LN-1, V( VPOS ),
+ CALL ZLARFX( 'Left', LM, LN-1, V( VPOS ),
$ DCONJG( TAU( TAUPOS ) ),
$ A( DPOS+NB-1, ST+1 ), LDA-1, WORK)
-
+
ENDIF
ENDIF
ENDIF
@@ -374,4 +374,4 @@
*
* END OF ZHB2ST_KERNELS
*
- END
+ END
diff --git a/lapack-netlib/SRC/zhecon_3.f b/lapack-netlib/SRC/zhecon_3.f
index 8c3a9f32b..9d2a240b6 100644
--- a/lapack-netlib/SRC/zhecon_3.f
+++ b/lapack-netlib/SRC/zhecon_3.f
@@ -19,7 +19,7 @@
* ===========
*
* SUBROUTINE ZHECON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND,
-* WORK, IWORK, INFO )
+* WORK, INFO )
*
* .. Scalar Arguments ..
* CHARACTER UPLO
@@ -27,7 +27,7 @@
* DOUBLE PRECISION ANORM, RCOND
* ..
* .. Array Arguments ..
-* INTEGER IPIV( * ), IWORK( * )
+* INTEGER IPIV( * )
* COMPLEX*16 A( LDA, * ), E ( * ), WORK( * )
* ..
*
@@ -129,11 +129,6 @@
*> WORK is COMPLEX*16 array, dimension (2*N)
*> \endverbatim
*>
-*> \param[out] IWORK
-*> \verbatim
-*> IWORK is INTEGER array, dimension (N)
-*> \endverbatim
-*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
diff --git a/lapack-netlib/SRC/zheevr.f b/lapack-netlib/SRC/zheevr.f
index 810373c83..def2d1f9d 100644
--- a/lapack-netlib/SRC/zheevr.f
+++ b/lapack-netlib/SRC/zheevr.f
@@ -210,7 +210,7 @@
*> eigenvalues are computed to high relative accuracy when
*> possible in future releases. The current code does not
*> make any guarantees about high relative accuracy, but
-*> furutre releases will. See J. Barlow and J. Demmel,
+*> future releases will. See J. Barlow and J. Demmel,
*> "Computing Accurate Eigensystems of Scaled Diagonally
*> Dominant Matrices", LAPACK Working Note #7, for a discussion
*> of which matrices define their eigenvalues to high relative
diff --git a/lapack-netlib/SRC/zheevr_2stage.f b/lapack-netlib/SRC/zheevr_2stage.f
index ab7f3374e..fe4a72160 100644
--- a/lapack-netlib/SRC/zheevr_2stage.f
+++ b/lapack-netlib/SRC/zheevr_2stage.f
@@ -217,7 +217,7 @@
*> eigenvalues are computed to high relative accuracy when
*> possible in future releases. The current code does not
*> make any guarantees about high relative accuracy, but
-*> furutre releases will. See J. Barlow and J. Demmel,
+*> future releases will. See J. Barlow and J. Demmel,
*> "Computing Accurate Eigensystems of Scaled Diagonally
*> Dominant Matrices", LAPACK Working Note #7, for a discussion
*> of which matrices define their eigenvalues to high relative
diff --git a/lapack-netlib/SRC/zhegs2.f b/lapack-netlib/SRC/zhegs2.f
index 0bdc653b9..aec526353 100644
--- a/lapack-netlib/SRC/zhegs2.f
+++ b/lapack-netlib/SRC/zhegs2.f
@@ -97,6 +97,7 @@
*> B is COMPLEX*16 array, dimension (LDB,N)
*> The triangular factor from the Cholesky factorization of B,
*> as returned by ZPOTRF.
+*> B is modified by the routine but restored on exit.
*> \endverbatim
*>
*> \param[in] LDB
diff --git a/lapack-netlib/SRC/zhegst.f b/lapack-netlib/SRC/zhegst.f
index d0c08a8f6..dcf5fe8b5 100644
--- a/lapack-netlib/SRC/zhegst.f
+++ b/lapack-netlib/SRC/zhegst.f
@@ -97,6 +97,7 @@
*> B is COMPLEX*16 array, dimension (LDB,N)
*> The triangular factor from the Cholesky factorization of B,
*> as returned by ZPOTRF.
+*> B is modified by the routine but restored on exit.
*> \endverbatim
*>
*> \param[in] LDB
diff --git a/lapack-netlib/SRC/zherfsx.f b/lapack-netlib/SRC/zherfsx.f
index d176b102c..fa11702a8 100644
--- a/lapack-netlib/SRC/zherfsx.f
+++ b/lapack-netlib/SRC/zherfsx.f
@@ -102,7 +102,7 @@
*> \param[in] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N)
-*> The symmetric matrix A. If UPLO = 'U', the leading N-by-N
+*> The Hermitian matrix A. If UPLO = 'U', the leading N-by-N
*> upper triangular part of A contains the upper triangular
*> part of the matrix A, and the strictly lower triangular
*> part of A is not referenced. If UPLO = 'L', the leading
@@ -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
@@ -306,14 +306,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.
@@ -321,9 +321,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/zhesv_aa.f b/lapack-netlib/SRC/zhesv_aa.f
index 8511f0e7d..5f1a9f4b3 100644
--- a/lapack-netlib/SRC/zhesv_aa.f
+++ b/lapack-netlib/SRC/zhesv_aa.f
@@ -42,7 +42,7 @@
*> matrices.
*>
*> Aasen's algorithm is used to factor A as
-*> A = U * T * U**H, if UPLO = 'U', or
+*> A = U**H * T * U, if UPLO = 'U', or
*> A = L * T * L**H, if UPLO = 'L',
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is Hermitian and tridiagonal. The factored form
@@ -86,7 +86,7 @@
*>
*> On exit, if INFO = 0, the tridiagonal matrix T and the
*> multipliers used to obtain the factor U or L from the
-*> factorization A = U*T*U**H or A = L*T*L**H as computed by
+*> factorization A = U**H*T*U or A = L*T*L**H as computed by
*> ZHETRF_AA.
*> \endverbatim
*>
@@ -230,7 +230,7 @@
RETURN
END IF
*
-* Compute the factorization A = U*T*U**H or A = L*T*L**H.
+* Compute the factorization A = U**H*T*U or A = L*T*L**H.
*
CALL ZHETRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
IF( INFO.EQ.0 ) THEN
diff --git a/lapack-netlib/SRC/zhesv_aa_2stage.f b/lapack-netlib/SRC/zhesv_aa_2stage.f
index ed221dc69..7a4e35f45 100644
--- a/lapack-netlib/SRC/zhesv_aa_2stage.f
+++ b/lapack-netlib/SRC/zhesv_aa_2stage.f
@@ -44,7 +44,7 @@
*> matrices.
*>
*> Aasen's 2-stage algorithm is used to factor A as
-*> A = U * T * U**H, if UPLO = 'U', or
+*> A = U**H * T * U, if UPLO = 'U', or
*> A = L * T * L**H, if UPLO = 'L',
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is Hermitian and band. The matrix T is
@@ -211,9 +211,7 @@
*
* .. Local Scalars ..
LOGICAL UPPER, TQUERY, WQUERY
- INTEGER I, J, K, I1, I2, TD
- INTEGER LDTB, LWKOPT, NB, KB, NT, IINFO
- COMPLEX PIV
+ INTEGER LWKOPT
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -263,7 +261,7 @@
RETURN
END IF
*
-* Compute the factorization A = U*T*U**H or A = L*T*L**H.
+* Compute the factorization A = U**H*T*U or A = L*T*L**H.
*
CALL ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
$ WORK, LWORK, INFO )
diff --git a/lapack-netlib/SRC/zhesvxx.f b/lapack-netlib/SRC/zhesvxx.f
index 375fc072d..20168185c 100644
--- a/lapack-netlib/SRC/zhesvxx.f
+++ b/lapack-netlib/SRC/zhesvxx.f
@@ -46,7 +46,7 @@
*>
*> ZHESVXX uses the diagonal pivoting factorization to compute the
*> solution to a complex*16 system of linear equations A * X = B, where
-*> A is an N-by-N symmetric matrix and X and B are N-by-NRHS
+*> A is an N-by-N Hermitian matrix and X and B are N-by-NRHS
*> matrices.
*>
*> If requested, both normwise and maximum componentwise error bounds
@@ -88,7 +88,7 @@
*> A = L * D * L**T, if UPLO = 'L',
*>
*> where U (or L) is a product of permutation and unit upper (lower)
-*> triangular matrices, and D is symmetric and block diagonal with
+*> triangular matrices, and D is Hermitian and block diagonal with
*> 1-by-1 and 2-by-2 diagonal blocks.
*>
*> 3. If some D(i,i)=0, so that D is exactly singular, then the
@@ -161,7 +161,7 @@
*> \param[in,out] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N)
-*> The symmetric matrix A. If UPLO = 'U', the leading N-by-N
+*> The Hermitian matrix A. If UPLO = 'U', the leading N-by-N
*> upper triangular part of A contains the upper triangular
*> part of the matrix A, and the strictly lower triangular
*> part of A is not referenced. If UPLO = 'L', the leading
@@ -378,7 +378,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
@@ -414,14 +414,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.
@@ -429,9 +429,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/zhetf2_rk.f b/lapack-netlib/SRC/zhetf2_rk.f
index 84d3a0248..6578214df 100644
--- a/lapack-netlib/SRC/zhetf2_rk.f
+++ b/lapack-netlib/SRC/zhetf2_rk.f
@@ -322,7 +322,7 @@
*
* Factorize A as U*D*U**H using the upper triangle of A
*
-* Initilize the first entry of array E, where superdiagonal
+* Initialize the first entry of array E, where superdiagonal
* elements of D are stored
*
E( 1 ) = CZERO
@@ -676,7 +676,7 @@
*
* Factorize A as L*D*L**H using the lower triangle of A
*
-* Initilize the unused last entry of the subdiagonal array E.
+* Initialize the unused last entry of the subdiagonal array E.
*
E( N ) = CZERO
*
diff --git a/lapack-netlib/SRC/zhetrd_2stage.f b/lapack-netlib/SRC/zhetrd_2stage.f
index 9d6a426a3..1a2c00a2f 100644
--- a/lapack-netlib/SRC/zhetrd_2stage.f
+++ b/lapack-netlib/SRC/zhetrd_2stage.f
@@ -123,23 +123,22 @@
*>
*> \param[out] HOUS2
*> \verbatim
-*> HOUS2 is COMPLEX*16 array, dimension LHOUS2, that
-*> store the Householder representation of the stage2
+*> HOUS2 is COMPLEX*16 array, dimension (LHOUS2)
+*> Stores the Householder representation of the stage2
*> band to tridiagonal.
*> \endverbatim
*>
*> \param[in] LHOUS2
*> \verbatim
*> LHOUS2 is INTEGER
-*> The dimension of the array HOUS2. LHOUS2 = MAX(1, dimension)
-*> If LWORK = -1, or LHOUS2=-1,
+*> The dimension of the array HOUS2.
+*> If LWORK = -1, or LHOUS2 = -1,
*> then a query is assumed; the routine
*> only calculates the optimal size of the HOUS2 array, returns
*> this value as the first entry of the HOUS2 array, and no error
*> message related to LHOUS2 is issued by XERBLA.
-*> LHOUS2 = MAX(1, dimension) where
-*> dimension = 4*N if VECT='N'
-*> not available now if VECT='H'
+*> If VECT='N', LHOUS2 = max(1, 4*n);
+*> if VECT='V', option not yet available.
*> \endverbatim
*>
*> \param[out] WORK
diff --git a/lapack-netlib/SRC/zhetrd_hb2st.F b/lapack-netlib/SRC/zhetrd_hb2st.F
index 86122cccc..b988f25ba 100644
--- a/lapack-netlib/SRC/zhetrd_hb2st.F
+++ b/lapack-netlib/SRC/zhetrd_hb2st.F
@@ -50,9 +50,9 @@
* Arguments:
* ==========
*
-*> \param[in] STAGE
+*> \param[in] STAGE1
*> \verbatim
-*> STAGE is CHARACTER*1
+*> STAGE1 is CHARACTER*1
*> = 'N': "No": to mention that the stage 1 of the reduction
*> from dense to band using the zhetrd_he2hb routine
*> was not called before this routine to reproduce AB.
@@ -512,8 +512,7 @@ C END IF
*
* Call the kernel
*
-#if defined(_OPENMP) && _OPENMP >= 201307
-
+#if defined(_OPENMP)
IF( TTYPE.NE.1 ) THEN
!$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1))
!$OMP$ DEPEND(in:WORK(MYID-1))
diff --git a/lapack-netlib/SRC/zhetrd_he2hb.f b/lapack-netlib/SRC/zhetrd_he2hb.f
index e33bf4b2b..b85b3889a 100644
--- a/lapack-netlib/SRC/zhetrd_he2hb.f
+++ b/lapack-netlib/SRC/zhetrd_he2hb.f
@@ -363,7 +363,7 @@
*
*
* Set the workspace of the triangular matrix T to zero once such a
-* way everytime T is generated the upper/lower portion will be always zero
+* way every time T is generated the upper/lower portion will be always zero
*
CALL ZLASET( "A", LDT, KD, ZERO, ZERO, WORK( TPOS ), LDT )
*
diff --git a/lapack-netlib/SRC/zhetrf_aa.f b/lapack-netlib/SRC/zhetrf_aa.f
index e355aed14..b80a84118 100644
--- a/lapack-netlib/SRC/zhetrf_aa.f
+++ b/lapack-netlib/SRC/zhetrf_aa.f
@@ -37,7 +37,7 @@
*> ZHETRF_AA computes the factorization of a complex hermitian matrix A
*> using the Aasen's algorithm. The form of the factorization is
*>
-*> A = U*T*U**H or A = L*T*L**H
+*> A = U**H*T*U or A = L*T*L**H
*>
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a hermitian tridiagonal matrix.
@@ -223,7 +223,7 @@
IF( UPPER ) THEN
*
* .....................................................
-* Factorize A as L*D*L**H using the upper triangle of A
+* Factorize A as U**H*D*U using the upper triangle of A
* .....................................................
*
* copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N))
@@ -256,7 +256,7 @@
$ A( MAX(1, J), J+1 ), LDA,
$ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) )
*
-* Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot)
+* Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
*
DO J2 = J+2, MIN(N, J+JB+1)
IPIV( J2 ) = IPIV( J2 ) + J
@@ -376,7 +376,7 @@
$ A( J+1, MAX(1, J) ), LDA,
$ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) )
*
-* Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot)
+* Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
*
DO J2 = J+2, MIN(N, J+JB+1)
IPIV( J2 ) = IPIV( J2 ) + J
diff --git a/lapack-netlib/SRC/zhetrf_aa_2stage.f b/lapack-netlib/SRC/zhetrf_aa_2stage.f
index 73c0ebe9a..f63713664 100644
--- a/lapack-netlib/SRC/zhetrf_aa_2stage.f
+++ b/lapack-netlib/SRC/zhetrf_aa_2stage.f
@@ -38,7 +38,7 @@
*> ZHETRF_AA_2STAGE computes the factorization of a double hermitian matrix A
*> using the Aasen's algorithm. The form of the factorization is
*>
-*> A = U*T*U**T or A = L*T*L**T
+*> A = U**H*T*U or A = L*T*L**H
*>
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a hermitian band matrix with the
@@ -66,7 +66,7 @@
*>
*> \param[in,out] A
*> \verbatim
-*> A is COMPLEX array, dimension (LDA,N)
+*> A is COMPLEX*16 array, dimension (LDA,N)
*> On entry, the hermitian matrix A. If UPLO = 'U', the leading
*> N-by-N upper triangular part of A contains the upper
*> triangular part of the matrix A, and the strictly lower
@@ -87,7 +87,7 @@
*>
*> \param[out] TB
*> \verbatim
-*> TB is COMPLEX array, dimension (LTB)
+*> TB is COMPLEX*16 array, dimension (LTB)
*> On exit, details of the LU factorization of the band matrix.
*> \endverbatim
*>
@@ -121,7 +121,7 @@
*>
*> \param[out] WORK
*> \verbatim
-*> WORK is COMPLEX workspace of size LWORK
+*> WORK is COMPLEX*16 workspace of size LWORK
*> \endverbatim
*>
*> \param[in] LWORK
@@ -276,7 +276,7 @@
IF( UPPER ) THEN
*
* .....................................................
-* Factorize A as L*D*L**T using the upper triangle of A
+* Factorize A as U**H*D*U using the upper triangle of A
* .....................................................
*
DO J = 0, NT-1
@@ -452,14 +452,17 @@ c END IF
* > Apply pivots to previous columns of L
CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
$ A( (J+1)*NB+1, I2 ), 1 )
-* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
- CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
- $ A( I1+1, I2 ), 1 )
+* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
+ IF( I2.GT.(I1+1) ) THEN
+ CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
+ $ A( I1+1, I2 ), 1 )
+ CALL ZLACGV( I2-I1-1, A( I1+1, I2 ), 1 )
+ END IF
CALL ZLACGV( I2-I1, A( I1, I1+1 ), LDA )
- CALL ZLACGV( I2-I1-1, A( I1+1, I2 ), 1 )
* > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
- CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA,
- $ A( I2, I2+1 ), LDA )
+ IF( I2.LT.N )
+ $ CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA,
+ $ A( I2, I2+1 ), LDA )
* > Swap A(I1, I1) with A(I2, I2)
PIV = A( I1, I1 )
A( I1, I1 ) = A( I2, I2 )
@@ -476,7 +479,7 @@ c END IF
ELSE
*
* .....................................................
-* Factorize A as L*D*L**T using the lower triangle of A
+* Factorize A as L*D*L**H using the lower triangle of A
* .....................................................
*
DO J = 0, NT-1
@@ -629,14 +632,17 @@ c END IF
* > Apply pivots to previous columns of L
CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
$ A( I2, (J+1)*NB+1 ), LDA )
-* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
- CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1,
- $ A( I2, I1+1 ), LDA )
+* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
+ IF( I2.GT.(I1+1) ) THEN
+ CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1,
+ $ A( I2, I1+1 ), LDA )
+ CALL ZLACGV( I2-I1-1, A( I2, I1+1 ), LDA )
+ END IF
CALL ZLACGV( I2-I1, A( I1+1, I1 ), 1 )
- CALL ZLACGV( I2-I1-1, A( I2, I1+1 ), LDA )
* > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
- CALL ZSWAP( N-I2, A( I2+1, I1 ), 1,
- $ A( I2+1, I2 ), 1 )
+ IF( I2.LT.N )
+ $ CALL ZSWAP( N-I2, A( I2+1, I1 ), 1,
+ $ A( I2+1, I2 ), 1 )
* > Swap A(I1, I1) with A(I2, I2)
PIV = A( I1, I1 )
A( I1, I1 ) = A( I2, I2 )
diff --git a/lapack-netlib/SRC/zhetri2.f b/lapack-netlib/SRC/zhetri2.f
index a7acff49f..ae43b14fe 100644
--- a/lapack-netlib/SRC/zhetri2.f
+++ b/lapack-netlib/SRC/zhetri2.f
@@ -62,7 +62,7 @@
*> \param[in,out] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N)
-*> On entry, the NB diagonal matrix D and the multipliers
+*> On entry, the block diagonal matrix D and the multipliers
*> used to obtain the factor U or L as computed by ZHETRF.
*>
*> On exit, if INFO = 0, the (symmetric) inverse of the original
@@ -82,7 +82,7 @@
*> \param[in] IPIV
*> \verbatim
*> IPIV is INTEGER array, dimension (N)
-*> Details of the interchanges and the NB structure of D
+*> Details of the interchanges and the block structure of D
*> as determined by ZHETRF.
*> \endverbatim
*>
diff --git a/lapack-netlib/SRC/zhetrs_aa.f b/lapack-netlib/SRC/zhetrs_aa.f
index 9d302b9cd..4b3253abc 100644
--- a/lapack-netlib/SRC/zhetrs_aa.f
+++ b/lapack-netlib/SRC/zhetrs_aa.f
@@ -38,8 +38,8 @@
*> \verbatim
*>
*> ZHETRS_AA solves a system of linear equations A*X = B with a complex
-*> hermitian matrix A using the factorization A = U*T*U**H or
-*> A = L*T*L**T computed by ZHETRF_AA.
+*> hermitian matrix A using the factorization A = U**H*T*U or
+*> A = L*T*L**H computed by ZHETRF_AA.
*> \endverbatim
*
* Arguments:
@@ -50,7 +50,7 @@
*> UPLO is CHARACTER*1
*> Specifies whether the details of the factorization are stored
*> as an upper or lower triangular matrix.
-*> = 'U': Upper triangular, form is A = U*T*U**H;
+*> = 'U': Upper triangular, form is A = U**H*T*U;
*> = 'L': Lower triangular, form is A = L*T*L**H.
*> \endverbatim
*>
@@ -98,14 +98,16 @@
*> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
-*> WORK is DOUBLE array, dimension (MAX(1,LWORK))
+*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
-*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2).
+*> LWORK is INTEGER
+*> The dimension of the array WORK. LWORK >= max(1,3*N-2).
+*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
@@ -199,61 +201,80 @@
*
IF( UPPER ) THEN
*
-* Solve A*X = B, where A = U*T*U**T.
+* Solve A*X = B, where A = U**H*T*U.
*
-* Pivot, P**T * B
+* 1) Forward substitution with U**H
*
- DO K = 1, N
- KP = IPIV( K )
- IF( KP.NE.K )
- $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
- END DO
+ IF( N.GT.1 ) THEN
*
-* Compute (U \P**T * B) -> B [ (U \P**T * B) ]
+* Pivot, P**T * B -> B
*
- CALL ZTRSM('L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
- $ B( 2, 1 ), LDB)
+ DO K = 1, N
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END DO
*
-* Compute T \ B -> B [ T \ (U \P**T * B) ]
+* Compute U**H \ B -> B [ (U**H \P**T * B) ]
*
- CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1)
+ CALL ZTRSM( 'L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ),
+ $ LDA, B( 2, 1 ), LDB )
+ END IF
+*
+* 2) Solve with triangular matrix T
+*
+* Compute T \ B -> B [ T \ (U**H \P**T * B) ]
+*
+ CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1 )
IF( N.GT.1 ) THEN
CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 2*N ), 1)
- CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1)
+ CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1 )
CALL ZLACGV( N-1, WORK( 1 ), 1 )
END IF
- CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
- $ INFO)
+ CALL ZGTSV( N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
+ $ INFO )
*
-* Compute (U**T \ B) -> B [ U**T \ (T \ (U \P**T * B) ) ]
+* 3) Backward substitution with U
*
- CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
- $ B(2, 1), LDB)
+ IF( N.GT.1 ) THEN
*
-* Pivot, P * B [ P * (U**T \ (T \ (U \P**T * B) )) ]
+* Compute U \ B -> B [ U \ (T \ (U**H \P**T * B) ) ]
*
- DO K = N, 1, -1
- KP = IPIV( K )
- IF( KP.NE.K )
- $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
- END DO
+ CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
+ $ LDA, B(2, 1), LDB)
+*
+* Pivot, P * B [ P * (U**H \ (T \ (U \P**T * B) )) ]
+*
+ DO K = N, 1, -1
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END DO
+ END IF
*
ELSE
*
-* Solve A*X = B, where A = L*T*L**T.
+* Solve A*X = B, where A = L*T*L**H.
*
-* Pivot, P**T * B
+* 1) Forward substitution with L
*
- DO K = 1, N
- KP = IPIV( K )
- IF( KP.NE.K )
- $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
- END DO
+ IF( N.GT.1 ) THEN
*
-* Compute (L \P**T * B) -> B [ (L \P**T * B) ]
+* Pivot, P**T * B -> B
*
- CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
- $ B(2, 1), LDB)
+ DO K = 1, N
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END DO
+*
+* Compute L \ B -> B [ (L \P**T * B) ]
+*
+ CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ),
+ $ LDA, B(2, 1), LDB)
+ END IF
+*
+* 2) Solve with triangular matrix T
*
* Compute T \ B -> B [ T \ (L \P**T * B) ]
*
@@ -266,18 +287,23 @@
CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
$ INFO)
*
-* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ]
+* 3) Backward substitution with L**H
*
- CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
- $ B( 2, 1 ), LDB)
+ IF( N.GT.1 ) THEN
*
-* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ]
+* Compute L**H \ B -> B [ L**H \ (T \ (L \P**T * B) ) ]
*
- DO K = N, 1, -1
- KP = IPIV( K )
- IF( KP.NE.K )
- $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
- END DO
+ CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ),
+ $ LDA, B( 2, 1 ), LDB)
+*
+* Pivot, P * B [ P * (L**H \ (T \ (L \P**T * B) )) ]
+*
+ DO K = N, 1, -1
+ KP = IPIV( K )
+ IF( KP.NE.K )
+ $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END DO
+ END IF
*
END IF
*
diff --git a/lapack-netlib/SRC/zhetrs_aa_2stage.f b/lapack-netlib/SRC/zhetrs_aa_2stage.f
index 7fcee1118..c621bd571 100644
--- a/lapack-netlib/SRC/zhetrs_aa_2stage.f
+++ b/lapack-netlib/SRC/zhetrs_aa_2stage.f
@@ -38,8 +38,8 @@
*> \verbatim
*>
*> ZHETRS_AA_2STAGE solves a system of linear equations A*X = B with a
-*> hermitian matrix A using the factorization A = U*T*U**T or
-*> A = L*T*L**T computed by ZHETRF_AA_2STAGE.
+*> hermitian matrix A using the factorization A = U**H*T*U or
+*> A = L*T*L**H computed by ZHETRF_AA_2STAGE.
*> \endverbatim
*
* Arguments:
@@ -50,8 +50,8 @@
*> UPLO is CHARACTER*1
*> Specifies whether the details of the factorization are stored
*> as an upper or lower triangular matrix.
-*> = 'U': Upper triangular, form is A = U*T*U**T;
-*> = 'L': Lower triangular, form is A = L*T*L**T.
+*> = 'U': Upper triangular, form is A = U**H*T*U;
+*> = 'L': Lower triangular, form is A = L*T*L**H.
*> \endverbatim
*>
*> \param[in] N
@@ -210,33 +210,33 @@
*
IF( UPPER ) THEN
*
-* Solve A*X = B, where A = U*T*U**T.
+* Solve A*X = B, where A = U**H*T*U.
*
IF( N.GT.NB ) THEN
*
-* Pivot, P**T * B
+* Pivot, P**T * B -> B
*
CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
*
-* Compute (U**T \P**T * B) -> B [ (U**T \P**T * B) ]
+* Compute (U**H \ B) -> B [ (U**H \P**T * B) ]
*
CALL ZTRSM( 'L', 'U', 'C', 'U', N-NB, NRHS, ONE, A(1, NB+1),
$ LDA, B(NB+1, 1), LDB)
*
END IF
*
-* Compute T \ B -> B [ T \ (U**T \P**T * B) ]
+* Compute T \ B -> B [ T \ (U**H \P**T * B) ]
*
CALL ZGBTRS( 'N', N, NB, NB, NRHS, TB, LDTB, IPIV2, B, LDB,
$ INFO)
IF( N.GT.NB ) THEN
*
-* Compute (U \ B) -> B [ U \ (T \ (U**T \P**T * B) ) ]
+* Compute (U \ B) -> B [ U \ (T \ (U**H \P**T * B) ) ]
*
CALL ZTRSM( 'L', 'U', 'N', 'U', N-NB, NRHS, ONE, A(1, NB+1),
$ LDA, B(NB+1, 1), LDB)
*
-* Pivot, P * B [ P * (U \ (T \ (U**T \P**T * B) )) ]
+* Pivot, P * B -> B [ P * (U \ (T \ (U**H \P**T * B) )) ]
*
CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
*
@@ -244,15 +244,15 @@
*
ELSE
*
-* Solve A*X = B, where A = L*T*L**T.
+* Solve A*X = B, where A = L*T*L**H.
*
IF( N.GT.NB ) THEN
*
-* Pivot, P**T * B
+* Pivot, P**T * B -> B
*
CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
*
-* Compute (L \P**T * B) -> B [ (L \P**T * B) ]
+* Compute (L \ B) -> B [ (L \P**T * B) ]
*
CALL ZTRSM( 'L', 'L', 'N', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
$ LDA, B(NB+1, 1), LDB)
@@ -265,12 +265,12 @@
$ INFO)
IF( N.GT.NB ) THEN
*
-* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ]
+* Compute (L**H \ B) -> B [ L**H \ (T \ (L \P**T * B) ) ]
*
CALL ZTRSM( 'L', 'L', 'C', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
$ LDA, B(NB+1, 1), LDB)
*
-* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ]
+* Pivot, P * B -> B [ P * (L**H \ (T \ (L \P**T * B) )) ]
*
CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
*
diff --git a/lapack-netlib/SRC/zhseqr.f b/lapack-netlib/SRC/zhseqr.f
index 1e8134c39..2ee874dfd 100644
--- a/lapack-netlib/SRC/zhseqr.f
+++ b/lapack-netlib/SRC/zhseqr.f
@@ -69,7 +69,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
@@ -86,7 +86,7 @@
*> set by a previous call to ZGEBAL, and then passed to ZGEHRD
*> when the matrix output by ZGEBAL 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
*>
@@ -98,17 +98,17 @@
*> triangular matrix T from the Schur decomposition (the
*> Schur form). 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 below.)
+*> INFO > 0 is given under the description of INFO below.)
*>
*> Unlike earlier versions of ZHSEQR, 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] W
@@ -131,7 +131,7 @@
*> if INFO = 0, Z contains Q*Z.
*> Normally Q is the unitary matrix generated by ZUNGHR
*> after the call to ZGEHRD 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
*>
@@ -139,7 +139,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
@@ -152,7 +152,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
@@ -170,21 +170,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, ZHSEQR 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
+*> > 0: if INFO = i, ZHSEQR failed to compute all of
+*> the eigenvalues. Elements 1:ilo-1 and i+1:n of W
+*> 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)
*>
@@ -192,19 +192,19 @@
*> value of H is upper Hessenberg and 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 unitary 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 unitary 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
*
@@ -244,8 +244,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
@@ -323,8 +323,8 @@
PARAMETER ( NTINY = 11 )
*
* ==== NL allocates some local workspace to help small matrices
-* . through a rare ZLAHQR failure. NL .GT. NTINY = 11 is
-* . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
+* . through a rare ZLAHQR 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/zla_gbrcond_c.f b/lapack-netlib/SRC/zla_gbrcond_c.f
index 20109124b..5b2dc46fc 100644
--- a/lapack-netlib/SRC/zla_gbrcond_c.f
+++ b/lapack-netlib/SRC/zla_gbrcond_c.f
@@ -133,13 +133,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_gbrcond_x.f b/lapack-netlib/SRC/zla_gbrcond_x.f
index 7e6c12ea5..17e9eede7 100644
--- a/lapack-netlib/SRC/zla_gbrcond_x.f
+++ b/lapack-netlib/SRC/zla_gbrcond_x.f
@@ -126,13 +126,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_gbrfsx_extended.f b/lapack-netlib/SRC/zla_gbrfsx_extended.f
index 7a850f1aa..a22b5592e 100644
--- a/lapack-netlib/SRC/zla_gbrfsx_extended.f
+++ b/lapack-netlib/SRC/zla_gbrfsx_extended.f
@@ -65,19 +65,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
@@ -269,7 +269,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/zla_gercond_c.f b/lapack-netlib/SRC/zla_gercond_c.f
index e629f90e8..a1c0df588 100644
--- a/lapack-netlib/SRC/zla_gercond_c.f
+++ b/lapack-netlib/SRC/zla_gercond_c.f
@@ -22,7 +22,7 @@
* LDAF, IPIV, C, CAPPLY,
* INFO, WORK, RWORK )
*
-* .. Scalar Aguments ..
+* .. Scalar Arguments ..
* CHARACTER TRANS
* LOGICAL CAPPLY
* INTEGER N, LDA, LDAF, INFO
@@ -114,13 +114,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
@@ -148,7 +148,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
-* .. Scalar Aguments ..
+* .. Scalar Arguments ..
CHARACTER TRANS
LOGICAL CAPPLY
INTEGER N, LDA, LDAF, INFO
diff --git a/lapack-netlib/SRC/zla_gercond_x.f b/lapack-netlib/SRC/zla_gercond_x.f
index 244bf58a3..3aa63ea84 100644
--- a/lapack-netlib/SRC/zla_gercond_x.f
+++ b/lapack-netlib/SRC/zla_gercond_x.f
@@ -107,13 +107,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_gerfsx_extended.f b/lapack-netlib/SRC/zla_gerfsx_extended.f
index 2e93e265e..e42ffa8e2 100644
--- a/lapack-netlib/SRC/zla_gerfsx_extended.f
+++ b/lapack-netlib/SRC/zla_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/zla_hercond_c.f b/lapack-netlib/SRC/zla_hercond_c.f
index 61cfe95f1..7c933cc3c 100644
--- a/lapack-netlib/SRC/zla_hercond_c.f
+++ b/lapack-netlib/SRC/zla_hercond_c.f
@@ -111,13 +111,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_hercond_x.f b/lapack-netlib/SRC/zla_hercond_x.f
index 9c19b487d..ee283c0b5 100644
--- a/lapack-netlib/SRC/zla_hercond_x.f
+++ b/lapack-netlib/SRC/zla_hercond_x.f
@@ -104,13 +104,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_herfsx_extended.f b/lapack-netlib/SRC/zla_herfsx_extended.f
index 5b43a58b9..8329080ef 100644
--- a/lapack-netlib/SRC/zla_herfsx_extended.f
+++ b/lapack-netlib/SRC/zla_herfsx_extended.f
@@ -66,11 +66,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
@@ -254,7 +254,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/zla_herpvgrw.f b/lapack-netlib/SRC/zla_herpvgrw.f
index 557d6e830..d414c371f 100644
--- a/lapack-netlib/SRC/zla_herpvgrw.f
+++ b/lapack-netlib/SRC/zla_herpvgrw.f
@@ -102,7 +102,7 @@
*> as determined by ZHETRF.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (2*N)
*> \endverbatim
diff --git a/lapack-netlib/SRC/zla_porcond_c.f b/lapack-netlib/SRC/zla_porcond_c.f
index a74295b41..2e591dd09 100644
--- a/lapack-netlib/SRC/zla_porcond_c.f
+++ b/lapack-netlib/SRC/zla_porcond_c.f
@@ -103,13 +103,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_porcond_x.f b/lapack-netlib/SRC/zla_porcond_x.f
index 0b2c84f42..4f409544f 100644
--- a/lapack-netlib/SRC/zla_porcond_x.f
+++ b/lapack-netlib/SRC/zla_porcond_x.f
@@ -96,13 +96,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_porfsx_extended.f b/lapack-netlib/SRC/zla_porfsx_extended.f
index 85dd42780..169a9a5d4 100644
--- a/lapack-netlib/SRC/zla_porfsx_extended.f
+++ b/lapack-netlib/SRC/zla_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/zla_porpvgrw.f b/lapack-netlib/SRC/zla_porpvgrw.f
index cd71635ec..f669b2864 100644
--- a/lapack-netlib/SRC/zla_porpvgrw.f
+++ b/lapack-netlib/SRC/zla_porpvgrw.f
@@ -86,7 +86,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/zla_syrcond_c.f b/lapack-netlib/SRC/zla_syrcond_c.f
index be9d14bd0..ff44d6c3b 100644
--- a/lapack-netlib/SRC/zla_syrcond_c.f
+++ b/lapack-netlib/SRC/zla_syrcond_c.f
@@ -111,13 +111,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_syrcond_x.f b/lapack-netlib/SRC/zla_syrcond_x.f
index 2d0269092..53022bbfb 100644
--- a/lapack-netlib/SRC/zla_syrcond_x.f
+++ b/lapack-netlib/SRC/zla_syrcond_x.f
@@ -104,13 +104,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] RWORK
+*> \param[out] RWORK
*> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/zla_syrfsx_extended.f b/lapack-netlib/SRC/zla_syrfsx_extended.f
index a9716fd23..69844c94b 100644
--- a/lapack-netlib/SRC/zla_syrfsx_extended.f
+++ b/lapack-netlib/SRC/zla_syrfsx_extended.f
@@ -66,11 +66,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
@@ -254,7 +254,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/zla_syrpvgrw.f b/lapack-netlib/SRC/zla_syrpvgrw.f
index ccf4fc2d6..82c9f52f8 100644
--- a/lapack-netlib/SRC/zla_syrpvgrw.f
+++ b/lapack-netlib/SRC/zla_syrpvgrw.f
@@ -102,7 +102,7 @@
*> as determined by ZSYTRF.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (2*N)
*> \endverbatim
diff --git a/lapack-netlib/SRC/zla_wwaddw.f b/lapack-netlib/SRC/zla_wwaddw.f
index b4f9df332..f06113a95 100644
--- a/lapack-netlib/SRC/zla_wwaddw.f
+++ b/lapack-netlib/SRC/zla_wwaddw.f
@@ -36,7 +36,7 @@
*> ZLA_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:
diff --git a/lapack-netlib/SRC/zlahef_aa.f b/lapack-netlib/SRC/zlahef_aa.f
index 8bad4aba9..ddd1e9493 100644
--- a/lapack-netlib/SRC/zlahef_aa.f
+++ b/lapack-netlib/SRC/zlahef_aa.f
@@ -288,8 +288,9 @@
*
* Swap A(I1, I2+1:N) with A(I2, I2+1:N)
*
- CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
- $ A( J1+I2-1, I2+1 ), LDA )
+ IF( I2.LT.M )
+ $ CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
+ $ A( J1+I2-1, I2+1 ), LDA )
*
* Swap A(I1, I1) with A(I2,I2)
*
@@ -329,13 +330,15 @@
* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
*
- IF( A( K, J+1 ).NE.ZERO ) THEN
- ALPHA = ONE / A( K, J+1 )
- CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
- CALL ZSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
- ELSE
- CALL ZLASET( 'Full', 1, M-J-1, ZERO, ZERO,
- $ A( K, J+2 ), LDA)
+ IF( J.LT.(M-1) ) THEN
+ IF( A( K, J+1 ).NE.ZERO ) THEN
+ ALPHA = ONE / A( K, J+1 )
+ CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
+ CALL ZSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
+ ELSE
+ CALL ZLASET( 'Full', 1, M-J-1, ZERO, ZERO,
+ $ A( K, J+2 ), LDA)
+ END IF
END IF
END IF
J = J + 1
@@ -440,8 +443,9 @@
*
* Swap A(I2+1:N, I1) with A(I2+1:N, I2)
*
- CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
- $ A( I2+1, J1+I2-1 ), 1 )
+ IF( I2.LT.M )
+ $ CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
+ $ A( I2+1, J1+I2-1 ), 1 )
*
* Swap A(I1, I1) with A(I2, I2)
*
@@ -481,13 +485,15 @@
* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
*
- IF( A( J+1, K ).NE.ZERO ) THEN
- ALPHA = ONE / A( J+1, K )
- CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
- CALL ZSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
- ELSE
- CALL ZLASET( 'Full', M-J-1, 1, ZERO, ZERO,
- $ A( J+2, K ), LDA )
+ IF( J.LT.(M-1) ) THEN
+ IF( A( J+1, K ).NE.ZERO ) THEN
+ ALPHA = ONE / A( J+1, K )
+ CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
+ CALL ZSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
+ ELSE
+ CALL ZLASET( 'Full', M-J-1, 1, ZERO, ZERO,
+ $ A( J+2, K ), LDA )
+ END IF
END IF
END IF
J = J + 1
diff --git a/lapack-netlib/SRC/zlahef_rk.f b/lapack-netlib/SRC/zlahef_rk.f
index d8d54f4ce..6a8549cf5 100644
--- a/lapack-netlib/SRC/zlahef_rk.f
+++ b/lapack-netlib/SRC/zlahef_rk.f
@@ -331,7 +331,7 @@
* Factorize the trailing columns of A using the upper triangle
* of A and working backwards, and compute the matrix W = U12*D
* for use in updating A11 (note that conjg(W) is actually stored)
-* Initilize the first entry of array E, where superdiagonal
+* Initialize the first entry of array E, where superdiagonal
* elements of D are stored
*
E( 1 ) = CZERO
@@ -789,7 +789,7 @@
* of A and working forwards, and compute the matrix W = L21*D
* for use in updating A22 (note that conjg(W) is actually stored)
*
-* Initilize the unused last entry of the subdiagonal array E.
+* Initialize the unused last entry of the subdiagonal array E.
*
E( N ) = CZERO
*
diff --git a/lapack-netlib/SRC/zlahqr.f b/lapack-netlib/SRC/zlahqr.f
index 19015b3fa..0a8318874 100644
--- a/lapack-netlib/SRC/zlahqr.f
+++ b/lapack-netlib/SRC/zlahqr.f
@@ -138,26 +138,26 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0: successful exit
-*> .GT. 0: if INFO = i, ZLAHQR failed to compute all the
+*> = 0: successful exit
+*> > 0: if INFO = i, ZLAHQR failed to compute all the
*> eigenvalues ILO to IHI in a total of 30 iterations
*> per eigenvalue; elements i+1:ihi of W contain
*> those eigenvalues which have been successfully
*> computed.
*>
-*> If INFO .GT. 0 and WANTT is .FALSE., then on exit,
+*> If INFO > 0 and WANTT is .FALSE., then on exit,
*> the remaining unconverged eigenvalues are the
*> eigenvalues of the upper Hessenberg matrix
-*> rows and columns ILO thorugh INFO of the final,
+*> rows and columns ILO through INFO of the final,
*> output value of H.
*>
-*> If INFO .GT. 0 and WANTT is .TRUE., then on exit
+*> If INFO > 0 and WANTT is .TRUE., then on exit
*> (*) (initial value of H)*U = U*(final value of H)
-*> where U is an orthognal matrix. The final
+*> where U is an orthogonal matrix. The final
*> value of H is upper Hessenberg and triangular in
*> rows and columns INFO+1 through IHI.
*>
-*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit
+*> If INFO > 0 and WANTZ is .TRUE., then on exit
*> (final value of Z) = (initial value of Z)*U
*> where U is the orthogonal matrix in (*)
*> (regardless of the value of WANTT.)
diff --git a/lapack-netlib/SRC/zlamswlq.f b/lapack-netlib/SRC/zlamswlq.f
index 0e0b0a1da..f32f5667c 100644
--- a/lapack-netlib/SRC/zlamswlq.f
+++ b/lapack-netlib/SRC/zlamswlq.f
@@ -1,3 +1,4 @@
+*> \brief \b ZLAMSWLQ
*
* Definition:
* ===========
diff --git a/lapack-netlib/SRC/zlamtsqr.f b/lapack-netlib/SRC/zlamtsqr.f
index 1ee732425..034c45505 100644
--- a/lapack-netlib/SRC/zlamtsqr.f
+++ b/lapack-netlib/SRC/zlamtsqr.f
@@ -1,3 +1,4 @@
+*> \brief \b ZLAMTSQR
*
* Definition:
* ===========
diff --git a/lapack-netlib/SRC/zlangb.f b/lapack-netlib/SRC/zlangb.f
index 949bb2c01..e40a470fd 100644
--- a/lapack-netlib/SRC/zlangb.f
+++ b/lapack-netlib/SRC/zlangb.f
@@ -130,6 +130,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM
INTEGER KL, KU, LDAB, N
@@ -147,14 +148,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J, K, L
- DOUBLE PRECISION SCALE, SUM, VALUE, TEMP
+ DOUBLE PRECISION SUM, VALUE, TEMP
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -207,15 +211,22 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
*
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 90 J = 1, N
L = MAX( 1, J-KU )
K = KU + 1 - J + L
- CALL ZLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANGB = VALUE
diff --git a/lapack-netlib/SRC/zlange.f b/lapack-netlib/SRC/zlange.f
index 5407decef..8162786fb 100644
--- a/lapack-netlib/SRC/zlange.f
+++ b/lapack-netlib/SRC/zlange.f
@@ -120,6 +120,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM
INTEGER LDA, M, N
@@ -137,14 +138,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J
- DOUBLE PRECISION SCALE, SUM, VALUE, TEMP
+ DOUBLE PRECISION SUM, VALUE, TEMP
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT
@@ -196,13 +200,19 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
*
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 90 J = 1, N
- CALL ZLASSQ( M, A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( M, A( 1, J ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANGE = VALUE
diff --git a/lapack-netlib/SRC/zlanhb.f b/lapack-netlib/SRC/zlanhb.f
index b3717804f..16b5c117c 100644
--- a/lapack-netlib/SRC/zlanhb.f
+++ b/lapack-netlib/SRC/zlanhb.f
@@ -137,6 +137,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM, UPLO
INTEGER K, LDAB, N
@@ -154,14 +155,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J, L
- DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
+ DOUBLE PRECISION ABSA, SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN, SQRT
@@ -233,39 +237,57 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
+*
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
+*
+* Sum off-diagonals
*
- SCALE = ZERO
- SUM = ONE
IF( K.GT.0 ) THEN
IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ),
- $ 1, SCALE, SUM )
+ $ 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE
L = K + 1
ELSE
DO 120 J = 1, N - 1
- CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE
L = 1
END IF
- SUM = 2*SUM
+ SSQ( 2 ) = 2*SSQ( 2 )
ELSE
L = 1
END IF
+*
+* Sum diagonal
+*
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
DO 130 J = 1, N
IF( DBLE( AB( L, J ) ).NE.ZERO ) THEN
ABSA = ABS( DBLE( AB( L, J ) ) )
- IF( SCALE.LT.ABSA ) THEN
- SUM = ONE + SUM*( SCALE / ABSA )**2
- SCALE = ABSA
+ IF( COLSSQ( 1 ).LT.ABSA ) THEN
+ COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
+ COLSSQ( 1 ) = ABSA
ELSE
- SUM = SUM + ( ABSA / SCALE )**2
+ COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
END IF
END IF
130 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANHB = VALUE
diff --git a/lapack-netlib/SRC/zlanhe.f b/lapack-netlib/SRC/zlanhe.f
index 7c7f7f3be..5aef9a756 100644
--- a/lapack-netlib/SRC/zlanhe.f
+++ b/lapack-netlib/SRC/zlanhe.f
@@ -129,6 +129,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM, UPLO
INTEGER LDA, N
@@ -146,14 +147,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J
- DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
+ DOUBLE PRECISION ABSA, SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, SQRT
@@ -223,31 +227,48 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
+*
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
+*
+* Sum off-diagonals
*
- SCALE = ZERO
- SUM = ONE
IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N
- CALL ZLASSQ( J-1, A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( J-1, A( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE
ELSE
DO 120 J = 1, N - 1
- CALL ZLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( N-J, A( J+1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE
END IF
- SUM = 2*SUM
+ SSQ( 2 ) = 2*SSQ( 2 )
+*
+* Sum diagonal
+*
DO 130 I = 1, N
IF( DBLE( A( I, I ) ).NE.ZERO ) THEN
ABSA = ABS( DBLE( A( I, I ) ) )
- IF( SCALE.LT.ABSA ) THEN
- SUM = ONE + SUM*( SCALE / ABSA )**2
- SCALE = ABSA
+ IF( SSQ( 1 ).LT.ABSA ) THEN
+ SSQ( 2 ) = ONE + SSQ( 2 )*( SSQ( 1 ) / ABSA )**2
+ SSQ( 1 ) = ABSA
ELSE
- SUM = SUM + ( ABSA / SCALE )**2
+ SSQ( 2 ) = SSQ( 2 ) + ( ABSA / SSQ( 1 ) )**2
END IF
END IF
130 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANHE = VALUE
diff --git a/lapack-netlib/SRC/zlanhp.f b/lapack-netlib/SRC/zlanhp.f
index 9ded60746..d795aeca9 100644
--- a/lapack-netlib/SRC/zlanhp.f
+++ b/lapack-netlib/SRC/zlanhp.f
@@ -122,6 +122,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM, UPLO
INTEGER N
@@ -139,14 +140,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J, K
- DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
+ DOUBLE PRECISION ABSA, SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, SQRT
@@ -225,31 +229,48 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
+*
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
+*
+* Sum off-diagonals
*
- SCALE = ZERO
- SUM = ONE
K = 2
IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N
- CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + J
110 CONTINUE
ELSE
DO 120 J = 1, N - 1
- CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + N - J + 1
120 CONTINUE
END IF
- SUM = 2*SUM
+ SSQ( 2 ) = 2*SSQ( 2 )
+*
+* Sum diagonal
+*
K = 1
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
DO 130 I = 1, N
IF( DBLE( AP( K ) ).NE.ZERO ) THEN
ABSA = ABS( DBLE( AP( K ) ) )
- IF( SCALE.LT.ABSA ) THEN
- SUM = ONE + SUM*( SCALE / ABSA )**2
- SCALE = ABSA
+ IF( COLSSQ( 1 ).LT.ABSA ) THEN
+ COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
+ COLSSQ( 1 ) = ABSA
ELSE
- SUM = SUM + ( ABSA / SCALE )**2
+ COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
END IF
END IF
IF( LSAME( UPLO, 'U' ) ) THEN
@@ -258,7 +279,8 @@
K = K + N - I + 1
END IF
130 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANHP = VALUE
diff --git a/lapack-netlib/SRC/zlanhs.f b/lapack-netlib/SRC/zlanhs.f
index f2d36b304..bd8e86be9 100644
--- a/lapack-netlib/SRC/zlanhs.f
+++ b/lapack-netlib/SRC/zlanhs.f
@@ -114,6 +114,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM
INTEGER LDA, N
@@ -131,14 +132,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J
- DOUBLE PRECISION SCALE, SUM, VALUE
+ DOUBLE PRECISION SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT
@@ -190,13 +194,20 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
*
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 90 J = 1, N
- CALL ZLASSQ( MIN( N, J+1 ), A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( MIN( N, J+1 ), A( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANHS = VALUE
diff --git a/lapack-netlib/SRC/zlansb.f b/lapack-netlib/SRC/zlansb.f
index 3468c49b3..245dcaf4b 100644
--- a/lapack-netlib/SRC/zlansb.f
+++ b/lapack-netlib/SRC/zlansb.f
@@ -135,6 +135,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM, UPLO
INTEGER K, LDAB, N
@@ -152,14 +153,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J, L
- DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
+ DOUBLE PRECISION ABSA, SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -227,29 +231,47 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
+*
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
+*
+* Sum off-diagonals
*
- SCALE = ZERO
- SUM = ONE
IF( K.GT.0 ) THEN
IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ),
- $ 1, SCALE, SUM )
+ $ 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE
L = K + 1
ELSE
DO 120 J = 1, N - 1
- CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE
L = 1
END IF
- SUM = 2*SUM
+ SSQ( 2 ) = 2*SSQ( 2 )
ELSE
L = 1
END IF
- CALL ZLASSQ( N, AB( L, 1 ), LDAB, SCALE, SUM )
- VALUE = SCALE*SQRT( SUM )
+*
+* Sum diagonal
+*
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( N, AB( L, 1 ), LDAB, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANSB = VALUE
diff --git a/lapack-netlib/SRC/zlansp.f b/lapack-netlib/SRC/zlansp.f
index 84fb972bb..fa9220487 100644
--- a/lapack-netlib/SRC/zlansp.f
+++ b/lapack-netlib/SRC/zlansp.f
@@ -120,6 +120,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM, UPLO
INTEGER N
@@ -137,14 +138,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J, K
- DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
+ DOUBLE PRECISION ABSA, SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, DIMAG, SQRT
@@ -219,40 +223,57 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
+*
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
+*
+* Sum off-diagonals
*
- SCALE = ZERO
- SUM = ONE
K = 2
IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N
- CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + J
110 CONTINUE
ELSE
DO 120 J = 1, N - 1
- CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + N - J + 1
120 CONTINUE
END IF
- SUM = 2*SUM
+ SSQ( 2 ) = 2*SSQ( 2 )
+*
+* Sum diagonal
+*
K = 1
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
DO 130 I = 1, N
IF( DBLE( AP( K ) ).NE.ZERO ) THEN
ABSA = ABS( DBLE( AP( K ) ) )
- IF( SCALE.LT.ABSA ) THEN
- SUM = ONE + SUM*( SCALE / ABSA )**2
- SCALE = ABSA
+ IF( COLSSQ( 1 ).LT.ABSA ) THEN
+ COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
+ COLSSQ( 1 ) = ABSA
ELSE
- SUM = SUM + ( ABSA / SCALE )**2
+ COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
END IF
END IF
IF( DIMAG( AP( K ) ).NE.ZERO ) THEN
ABSA = ABS( DIMAG( AP( K ) ) )
- IF( SCALE.LT.ABSA ) THEN
- SUM = ONE + SUM*( SCALE / ABSA )**2
- SCALE = ABSA
+ IF( COLSSQ( 1 ).LT.ABSA ) THEN
+ COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
+ COLSSQ( 1 ) = ABSA
ELSE
- SUM = SUM + ( ABSA / SCALE )**2
+ COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
END IF
END IF
IF( LSAME( UPLO, 'U' ) ) THEN
@@ -261,7 +282,8 @@
K = K + N - I + 1
END IF
130 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANSP = VALUE
diff --git a/lapack-netlib/SRC/zlansy.f b/lapack-netlib/SRC/zlansy.f
index 58269a911..e022f85e1 100644
--- a/lapack-netlib/SRC/zlansy.f
+++ b/lapack-netlib/SRC/zlansy.f
@@ -128,6 +128,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM, UPLO
INTEGER LDA, N
@@ -145,14 +146,17 @@
* ..
* .. Local Scalars ..
INTEGER I, J
- DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
+ DOUBLE PRECISION ABSA, SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
@@ -218,21 +222,39 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
+*
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
+*
+* Sum off-diagonals
*
- SCALE = ZERO
- SUM = ONE
IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N
- CALL ZLASSQ( J-1, A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( J-1, A( 1, J ), 1, COLSSQ(1), COLSSQ(2) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE
ELSE
DO 120 J = 1, N - 1
- CALL ZLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( N-J, A( J+1, J ), 1, COLSSQ(1), COLSSQ(2) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE
END IF
- SUM = 2*SUM
- CALL ZLASSQ( N, A, LDA+1, SCALE, SUM )
- VALUE = SCALE*SQRT( SUM )
+ SSQ( 2 ) = 2*SSQ( 2 )
+*
+* Sum diagonal
+*
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( N, A, LDA+1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANSY = VALUE
diff --git a/lapack-netlib/SRC/zlantb.f b/lapack-netlib/SRC/zlantb.f
index 3077ba151..f02509223 100644
--- a/lapack-netlib/SRC/zlantb.f
+++ b/lapack-netlib/SRC/zlantb.f
@@ -146,6 +146,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER DIAG, NORM, UPLO
INTEGER K, LDAB, N
@@ -164,14 +165,17 @@
* .. Local Scalars ..
LOGICAL UDIAG
INTEGER I, J, L
- DOUBLE PRECISION SCALE, SUM, VALUE
+ DOUBLE PRECISION SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
@@ -313,46 +317,61 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
*
IF( LSAME( UPLO, 'U' ) ) THEN
IF( LSAME( DIAG, 'U' ) ) THEN
- SCALE = ONE
- SUM = N
+ SSQ( 1 ) = ONE
+ SSQ( 2 ) = N
IF( K.GT.0 ) THEN
DO 280 J = 2, N
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( J-1, K ),
- $ AB( MAX( K+2-J, 1 ), J ), 1, SCALE,
- $ SUM )
+ $ AB( MAX( K+2-J, 1 ), J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
280 CONTINUE
END IF
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 290 J = 1, N
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( J, K+1 ), AB( MAX( K+2-J, 1 ), J ),
- $ 1, SCALE, SUM )
+ $ 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
290 CONTINUE
END IF
ELSE
IF( LSAME( DIAG, 'U' ) ) THEN
- SCALE = ONE
- SUM = N
+ SSQ( 1 ) = ONE
+ SSQ( 2 ) = N
IF( K.GT.0 ) THEN
DO 300 J = 1, N - 1
- CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
300 CONTINUE
END IF
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 310 J = 1, N
- CALL ZLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
310 CONTINUE
END IF
END IF
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANTB = VALUE
diff --git a/lapack-netlib/SRC/zlantp.f b/lapack-netlib/SRC/zlantp.f
index 69dbaa5bc..d32a00f13 100644
--- a/lapack-netlib/SRC/zlantp.f
+++ b/lapack-netlib/SRC/zlantp.f
@@ -130,6 +130,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER DIAG, NORM, UPLO
INTEGER N
@@ -148,14 +149,17 @@
* .. Local Scalars ..
LOGICAL UDIAG
INTEGER I, J, K
- DOUBLE PRECISION SCALE, SUM, VALUE
+ DOUBLE PRECISION SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
@@ -308,45 +312,64 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
*
IF( LSAME( UPLO, 'U' ) ) THEN
IF( LSAME( DIAG, 'U' ) ) THEN
- SCALE = ONE
- SUM = N
+ SSQ( 1 ) = ONE
+ SSQ( 2 ) = N
K = 2
DO 280 J = 2, N
- CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( J-1, AP( K ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + J
280 CONTINUE
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
K = 1
DO 290 J = 1, N
- CALL ZLASSQ( J, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( J, AP( K ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + J
290 CONTINUE
END IF
ELSE
IF( LSAME( DIAG, 'U' ) ) THEN
- SCALE = ONE
- SUM = N
+ SSQ( 1 ) = ONE
+ SSQ( 2 ) = N
K = 2
DO 300 J = 1, N - 1
- CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( N-J, AP( K ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + N - J + 1
300 CONTINUE
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
K = 1
DO 310 J = 1, N
- CALL ZLASSQ( N-J+1, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( N-J+1, AP( K ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + N - J + 1
310 CONTINUE
END IF
END IF
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANTP = VALUE
diff --git a/lapack-netlib/SRC/zlantr.f b/lapack-netlib/SRC/zlantr.f
index 04ee482f7..7d63c972e 100644
--- a/lapack-netlib/SRC/zlantr.f
+++ b/lapack-netlib/SRC/zlantr.f
@@ -147,6 +147,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER DIAG, NORM, UPLO
INTEGER LDA, M, N
@@ -165,14 +166,17 @@
* .. Local Scalars ..
LOGICAL UDIAG
INTEGER I, J
- DOUBLE PRECISION SCALE, SUM, VALUE
+ DOUBLE PRECISION SUM, VALUE
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN
* ..
* .. External Subroutines ..
- EXTERNAL ZLASSQ
+ EXTERNAL ZLASSQ, DCOMBSSQ
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT
@@ -283,7 +287,7 @@
END IF
ELSE
IF( LSAME( DIAG, 'U' ) ) THEN
- DO 210 I = 1, N
+ DO 210 I = 1, MIN( M, N )
WORK( I ) = ONE
210 CONTINUE
DO 220 I = N + 1, M
@@ -313,38 +317,56 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
*
* Find normF(A).
+* SSQ(1) is scale
+* SSQ(2) is sum-of-squares
+* For better accuracy, sum each column separately.
*
IF( LSAME( UPLO, 'U' ) ) THEN
IF( LSAME( DIAG, 'U' ) ) THEN
- SCALE = ONE
- SUM = MIN( M, N )
+ SSQ( 1 ) = ONE
+ SSQ( 2 ) = MIN( M, N )
DO 290 J = 2, N
- CALL ZLASSQ( MIN( M, J-1 ), A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( MIN( M, J-1 ), A( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
290 CONTINUE
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 300 J = 1, N
- CALL ZLASSQ( MIN( M, J ), A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( MIN( M, J ), A( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
300 CONTINUE
END IF
ELSE
IF( LSAME( DIAG, 'U' ) ) THEN
- SCALE = ONE
- SUM = MIN( M, N )
+ SSQ( 1 ) = ONE
+ SSQ( 2 ) = MIN( M, N )
DO 310 J = 1, N
- CALL ZLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( M-J, A( MIN( M, J+1 ), J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
310 CONTINUE
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 320 J = 1, N
- CALL ZLASSQ( M-J+1, A( J, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL ZLASSQ( M-J+1, A( J, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL DCOMBSSQ( SSQ, COLSSQ )
320 CONTINUE
END IF
END IF
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
ZLANTR = VALUE
diff --git a/lapack-netlib/SRC/zlaqps.f b/lapack-netlib/SRC/zlaqps.f
index c142e8c69..66c721517 100644
--- a/lapack-netlib/SRC/zlaqps.f
+++ b/lapack-netlib/SRC/zlaqps.f
@@ -127,7 +127,7 @@
*> \param[in,out] AUXV
*> \verbatim
*> AUXV is COMPLEX*16 array, dimension (NB)
-*> Auxiliar vector.
+*> Auxiliary vector.
*> \endverbatim
*>
*> \param[in,out] F
diff --git a/lapack-netlib/SRC/zlaqr0.f b/lapack-netlib/SRC/zlaqr0.f
index 59b8ed7a6..feffe9782 100644
--- a/lapack-netlib/SRC/zlaqr0.f
+++ b/lapack-netlib/SRC/zlaqr0.f
@@ -66,7 +66,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
@@ -79,12 +79,12 @@
*> IHI is INTEGER
*>
*> It is assumed that H is already upper triangular in rows
-*> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
+*> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
*> previous call to ZGEBAL, and then passed to ZGEHRD when the
*> matrix output by ZGEBAL 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
*>
@@ -96,17 +96,17 @@
*> contains the upper triangular matrix T from the Schur
*> decomposition (the Schur form). If INFO = 0 and WANT is
*> .FALSE., then the contents of H are unspecified on exit.
-*> (The output value of H when INFO.GT.0 is given under the
+*> (The output value of H when INFO > 0 is given under the
*> description of INFO below.)
*>
-*> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
+*> This subroutine may explicitly set 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] W
@@ -128,7 +128,7 @@
*> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE..
-*> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
+*> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*> \endverbatim
*>
*> \param[in,out] Z
@@ -138,7 +138,7 @@
*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
-*> (The output value of Z when INFO.GT.0 is given under
+*> (The output value of Z when INFO > 0 is given under
*> the description of INFO below.)
*> \endverbatim
*>
@@ -146,7 +146,7 @@
*> \verbatim
*> LDZ is INTEGER
*> The leading dimension of the array Z. if WANTZ is .TRUE.
-*> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.
+*> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
*> \endverbatim
*>
*> \param[out] WORK
@@ -159,7 +159,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, but LWORK typically as large as 6*N may
*> be required for optimal performance. A workspace query
*> to determine the optimal workspace size is recommended.
@@ -175,19 +175,19 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0: successful exit
-*> .GT. 0: if INFO = i, ZLAQR0 failed to compute all of
+*> = 0: successful exit
+*> > 0: if INFO = i, ZLAQR0 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 WANT is .FALSE., then on exit,
+*> If INFO > 0 and WANT is .FALSE., 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 WANTT is .TRUE., then on exit
+*> If INFO > 0 and WANTT is .TRUE., then on exit
*>
*> (*) (initial value of H)*U = U*(final value of H)
*>
@@ -195,7 +195,7 @@
*> value of H is upper Hessenberg and triangular in
*> rows and columns INFO+1 through IHI.
*>
-*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit
+*> If INFO > 0 and WANTZ is .TRUE., then on exit
*>
*> (final value of Z(ILO:IHI,ILOZ:IHIZ)
*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
@@ -203,7 +203,7 @@
*> where U is the unitary matrix in (*) (regard-
*> less of the value of WANTT.)
*>
-*> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
+*> If INFO > 0 and WANTZ is .FALSE., then Z is not
*> accessed.
*> \endverbatim
*
@@ -641,7 +641,7 @@
END IF
END IF
*
-* ==== Use up to NS of the the smallest magnatiude
+* ==== Use up to NS of the the smallest magnitude
* . shifts. If there aren't NS shifts available,
* . then use them all, possibly dropping one to
* . make the number of shifts even. ====
diff --git a/lapack-netlib/SRC/zlaqr1.f b/lapack-netlib/SRC/zlaqr1.f
index 34341cb10..fc2df3cb4 100644
--- a/lapack-netlib/SRC/zlaqr1.f
+++ b/lapack-netlib/SRC/zlaqr1.f
@@ -64,7 +64,7 @@
*> \verbatim
*> LDH is INTEGER
*> The leading dimension of H as declared in
-*> the calling procedure. LDH.GE.N
+*> the calling procedure. LDH >= N
*> \endverbatim
*>
*> \param[in] S1
diff --git a/lapack-netlib/SRC/zlaqr2.f b/lapack-netlib/SRC/zlaqr2.f
index e6e2ea48c..b5434e899 100644
--- a/lapack-netlib/SRC/zlaqr2.f
+++ b/lapack-netlib/SRC/zlaqr2.f
@@ -103,7 +103,7 @@
*> \param[in] NW
*> \verbatim
*> NW is INTEGER
-*> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
+*> Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
*> \endverbatim
*>
*> \param[in,out] H
@@ -121,7 +121,7 @@
*> \verbatim
*> LDH is INTEGER
*> Leading dimension of H just as declared in the calling
-*> subroutine. N .LE. LDH
+*> subroutine. N <= LDH
*> \endverbatim
*>
*> \param[in] ILOZ
@@ -133,7 +133,7 @@
*> \verbatim
*> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be
-*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
+*> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
*> \endverbatim
*>
*> \param[in,out] Z
@@ -149,7 +149,7 @@
*> \verbatim
*> LDZ is INTEGER
*> The leading dimension of Z just as declared in the
-*> calling subroutine. 1 .LE. LDZ.
+*> calling subroutine. 1 <= LDZ.
*> \endverbatim
*>
*> \param[out] NS
@@ -186,13 +186,13 @@
*> \verbatim
*> LDV is INTEGER
*> The leading dimension of V just as declared in the
-*> calling subroutine. NW .LE. LDV
+*> calling subroutine. NW <= LDV
*> \endverbatim
*>
*> \param[in] NH
*> \verbatim
*> NH is INTEGER
-*> The number of columns of T. NH.GE.NW.
+*> The number of columns of T. NH >= NW.
*> \endverbatim
*>
*> \param[out] T
@@ -204,14 +204,14 @@
*> \verbatim
*> LDT is INTEGER
*> The leading dimension of T just as declared in the
-*> calling subroutine. NW .LE. LDT
+*> calling subroutine. NW <= LDT
*> \endverbatim
*>
*> \param[in] NV
*> \verbatim
*> NV is INTEGER
*> The number of rows of work array WV available for
-*> workspace. NV.GE.NW.
+*> workspace. NV >= NW.
*> \endverbatim
*>
*> \param[out] WV
@@ -223,7 +223,7 @@
*> \verbatim
*> LDWV is INTEGER
*> The leading dimension of W just as declared in the
-*> calling subroutine. NW .LE. LDV
+*> calling subroutine. NW <= LDV
*> \endverbatim
*>
*> \param[out] WORK
diff --git a/lapack-netlib/SRC/zlaqr3.f b/lapack-netlib/SRC/zlaqr3.f
index 64ab59f31..dfb798ca9 100644
--- a/lapack-netlib/SRC/zlaqr3.f
+++ b/lapack-netlib/SRC/zlaqr3.f
@@ -100,7 +100,7 @@
*> \param[in] NW
*> \verbatim
*> NW is INTEGER
-*> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
+*> Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
*> \endverbatim
*>
*> \param[in,out] H
@@ -118,7 +118,7 @@
*> \verbatim
*> LDH is INTEGER
*> Leading dimension of H just as declared in the calling
-*> subroutine. N .LE. LDH
+*> subroutine. N <= LDH
*> \endverbatim
*>
*> \param[in] ILOZ
@@ -130,7 +130,7 @@
*> \verbatim
*> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be
-*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
+*> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
*> \endverbatim
*>
*> \param[in,out] Z
@@ -146,7 +146,7 @@
*> \verbatim
*> LDZ is INTEGER
*> The leading dimension of Z just as declared in the
-*> calling subroutine. 1 .LE. LDZ.
+*> calling subroutine. 1 <= LDZ.
*> \endverbatim
*>
*> \param[out] NS
@@ -183,13 +183,13 @@
*> \verbatim
*> LDV is INTEGER
*> The leading dimension of V just as declared in the
-*> calling subroutine. NW .LE. LDV
+*> calling subroutine. NW <= LDV
*> \endverbatim
*>
*> \param[in] NH
*> \verbatim
*> NH is INTEGER
-*> The number of columns of T. NH.GE.NW.
+*> The number of columns of T. NH >= NW.
*> \endverbatim
*>
*> \param[out] T
@@ -201,14 +201,14 @@
*> \verbatim
*> LDT is INTEGER
*> The leading dimension of T just as declared in the
-*> calling subroutine. NW .LE. LDT
+*> calling subroutine. NW <= LDT
*> \endverbatim
*>
*> \param[in] NV
*> \verbatim
*> NV is INTEGER
*> The number of rows of work array WV available for
-*> workspace. NV.GE.NW.
+*> workspace. NV >= NW.
*> \endverbatim
*>
*> \param[out] WV
@@ -220,7 +220,7 @@
*> \verbatim
*> LDWV is INTEGER
*> The leading dimension of W just as declared in the
-*> calling subroutine. NW .LE. LDV
+*> calling subroutine. NW <= LDV
*> \endverbatim
*>
*> \param[out] WORK
diff --git a/lapack-netlib/SRC/zlaqr4.f b/lapack-netlib/SRC/zlaqr4.f
index 012fa37e2..a88f6508e 100644
--- a/lapack-netlib/SRC/zlaqr4.f
+++ b/lapack-netlib/SRC/zlaqr4.f
@@ -73,7 +73,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
@@ -85,12 +85,12 @@
*> \verbatim
*> IHI is INTEGER
*> It is assumed that H is already upper triangular in rows
-*> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
+*> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
*> previous call to ZGEBAL, and then passed to ZGEHRD when the
*> matrix output by ZGEBAL 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
*>
@@ -102,17 +102,17 @@
*> contains the upper triangular matrix T from the Schur
*> decomposition (the Schur form). If INFO = 0 and WANT is
*> .FALSE., then the contents of H are unspecified on exit.
-*> (The output value of H when INFO.GT.0 is given under the
+*> (The output value of H when INFO > 0 is given under the
*> description of INFO below.)
*>
-*> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
+*> This subroutine may explicitly set 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] W
@@ -134,7 +134,7 @@
*> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE..
-*> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
+*> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*> \endverbatim
*>
*> \param[in,out] Z
@@ -144,7 +144,7 @@
*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
-*> (The output value of Z when INFO.GT.0 is given under
+*> (The output value of Z when INFO > 0 is given under
*> the description of INFO below.)
*> \endverbatim
*>
@@ -152,7 +152,7 @@
*> \verbatim
*> LDZ is INTEGER
*> The leading dimension of the array Z. if WANTZ is .TRUE.
-*> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.
+*> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
*> \endverbatim
*>
*> \param[out] WORK
@@ -165,7 +165,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, but LWORK typically as large as 6*N may
*> be required for optimal performance. A workspace query
*> to determine the optimal workspace size is recommended.
@@ -182,18 +182,18 @@
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
-*> .GT. 0: if INFO = i, ZLAQR4 failed to compute all of
+*> > 0: if INFO = i, ZLAQR4 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 WANT is .FALSE., then on exit,
+*> If INFO > 0 and WANT is .FALSE., 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 WANTT is .TRUE., then on exit
+*> If INFO > 0 and WANTT is .TRUE., then on exit
*>
*> (*) (initial value of H)*U = U*(final value of H)
*>
@@ -201,7 +201,7 @@
*> value of H is upper Hessenberg and triangular in
*> rows and columns INFO+1 through IHI.
*>
-*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit
+*> If INFO > 0 and WANTZ is .TRUE., then on exit
*>
*> (final value of Z(ILO:IHI,ILOZ:IHIZ)
*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
@@ -209,7 +209,7 @@
*> where U is the unitary matrix in (*) (regard-
*> less of the value of WANTT.)
*>
-*> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
+*> If INFO > 0 and WANTZ is .FALSE., then Z is not
*> accessed.
*> \endverbatim
*
@@ -641,7 +641,7 @@
END IF
END IF
*
-* ==== Use up to NS of the the smallest magnatiude
+* ==== Use up to NS of the the smallest magnitude
* . shifts. If there aren't NS shifts available,
* . then use them all, possibly dropping one to
* . make the number of shifts even. ====
diff --git a/lapack-netlib/SRC/zlaqr5.f b/lapack-netlib/SRC/zlaqr5.f
index 0dfbce82c..9ff7e7eca 100644
--- a/lapack-netlib/SRC/zlaqr5.f
+++ b/lapack-netlib/SRC/zlaqr5.f
@@ -125,7 +125,7 @@
*> \verbatim
*> LDH is INTEGER
*> LDH is the leading dimension of H just as declared in the
-*> calling procedure. LDH.GE.MAX(1,N).
+*> calling procedure. LDH >= MAX(1,N).
*> \endverbatim
*>
*> \param[in] ILOZ
@@ -137,7 +137,7 @@
*> \verbatim
*> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be
-*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
+*> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N
*> \endverbatim
*>
*> \param[in,out] Z
@@ -153,7 +153,7 @@
*> \verbatim
*> LDZ is INTEGER
*> LDA is the leading dimension of Z just as declared in
-*> the calling procedure. LDZ.GE.N.
+*> the calling procedure. LDZ >= N.
*> \endverbatim
*>
*> \param[out] V
@@ -165,7 +165,7 @@
*> \verbatim
*> LDV is INTEGER
*> LDV is the leading dimension of V as declared in the
-*> calling procedure. LDV.GE.3.
+*> calling procedure. LDV >= 3.
*> \endverbatim
*>
*> \param[out] U
@@ -177,33 +177,14 @@
*> \verbatim
*> LDU is INTEGER
*> LDU is the leading dimension of U just as declared in the
-*> in the calling subroutine. LDU.GE.3*NSHFTS-3.
-*> \endverbatim
-*>
-*> \param[in] NH
-*> \verbatim
-*> NH is INTEGER
-*> NH is the number of columns in array WH available for
-*> workspace. NH.GE.1.
-*> \endverbatim
-*>
-*> \param[out] WH
-*> \verbatim
-*> WH is COMPLEX*16 array, dimension (LDWH,NH)
-*> \endverbatim
-*>
-*> \param[in] LDWH
-*> \verbatim
-*> LDWH is INTEGER
-*> Leading dimension of WH just as declared in the
-*> calling procedure. LDWH.GE.3*NSHFTS-3.
+*> in the calling subroutine. LDU >= 3*NSHFTS-3.
*> \endverbatim
*>
*> \param[in] NV
*> \verbatim
*> NV is INTEGER
*> NV is the number of rows in WV agailable for workspace.
-*> NV.GE.1.
+*> NV >= 1.
*> \endverbatim
*>
*> \param[out] WV
@@ -215,9 +196,28 @@
*> \verbatim
*> LDWV is INTEGER
*> LDWV is the leading dimension of WV as declared in the
-*> in the calling subroutine. LDWV.GE.NV.
+*> in the calling subroutine. LDWV >= NV.
*> \endverbatim
*
+*> \param[in] NH
+*> \verbatim
+*> NH is INTEGER
+*> NH is the number of columns in array WH available for
+*> workspace. NH >= 1.
+*> \endverbatim
+*>
+*> \param[out] WH
+*> \verbatim
+*> WH is COMPLEX*16 array, dimension (LDWH,NH)
+*> \endverbatim
+*>
+*> \param[in] LDWH
+*> \verbatim
+*> LDWH is INTEGER
+*> Leading dimension of WH just as declared in the
+*> calling procedure. LDWH >= 3*NSHFTS-3.
+*> \endverbatim
+*>
* Authors:
* ========
*
diff --git a/lapack-netlib/SRC/zlarfb.f b/lapack-netlib/SRC/zlarfb.f
index b4a2b4d1a..3da49f2fc 100644
--- a/lapack-netlib/SRC/zlarfb.f
+++ b/lapack-netlib/SRC/zlarfb.f
@@ -92,6 +92,8 @@
*> K is INTEGER
*> The order of the matrix T (= the number of elementary
*> reflectors whose product defines the block reflector).
+*> If SIDE = 'L', M >= K >= 0;
+*> if SIDE = 'R', N >= K >= 0.
*> \endverbatim
*>
*> \param[in] V
diff --git a/lapack-netlib/SRC/zlarfx.f b/lapack-netlib/SRC/zlarfx.f
index 685d164eb..ba6d4ed74 100644
--- a/lapack-netlib/SRC/zlarfx.f
+++ b/lapack-netlib/SRC/zlarfx.f
@@ -94,7 +94,7 @@
*> \param[in] LDC
*> \verbatim
*> LDC is INTEGER
-*> The leading dimension of the array C. LDA >= max(1,M).
+*> The leading dimension of the array C. LDC >= max(1,M).
*> \endverbatim
*>
*> \param[out] WORK
diff --git a/lapack-netlib/SRC/zlarfy.f b/lapack-netlib/SRC/zlarfy.f
index 57605731b..4c9e08bac 100644
--- a/lapack-netlib/SRC/zlarfy.f
+++ b/lapack-netlib/SRC/zlarfy.f
@@ -103,7 +103,7 @@
*
*> \date December 2016
*
-*> \ingroup complex16_eig
+*> \ingroup complex16OTHERauxiliary
*
* =====================================================================
SUBROUTINE ZLARFY( UPLO, N, V, INCV, TAU, C, LDC, WORK )
diff --git a/lapack-netlib/SRC/zlarrv.f b/lapack-netlib/SRC/zlarrv.f
index 67a67584c..23976dbef 100644
--- a/lapack-netlib/SRC/zlarrv.f
+++ b/lapack-netlib/SRC/zlarrv.f
@@ -143,7 +143,7 @@
*> RTOL2 is DOUBLE PRECISION
*> Parameters for bisection.
*> An interval [LEFT,RIGHT] has converged if
-*> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
+*> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*> \endverbatim
*>
*> \param[in,out] W
diff --git a/lapack-netlib/SRC/zlassq.f b/lapack-netlib/SRC/zlassq.f
index fd13811bd..dccec988d 100644
--- a/lapack-netlib/SRC/zlassq.f
+++ b/lapack-netlib/SRC/zlassq.f
@@ -41,7 +41,7 @@
*> where x( i ) = abs( X( 1 + ( i - 1 )*INCX ) ). The value of sumsq is
*> assumed to be at least unity and the value of ssq will then satisfy
*>
-*> 1.0 .le. ssq .le. ( sumsq + 2*n ).
+*> 1.0 <= ssq <= ( sumsq + 2*n ).
*>
*> scale is assumed to be non-negative and scl returns the value
*>
@@ -65,7 +65,7 @@
*>
*> \param[in] X
*> \verbatim
-*> X is COMPLEX*16 array, dimension (N)
+*> X is COMPLEX*16 array, dimension (1+(N-1)*INCX)
*> The vector x as described above.
*> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
*> \endverbatim
diff --git a/lapack-netlib/SRC/zlaswlq.f b/lapack-netlib/SRC/zlaswlq.f
index 24dd41d79..990630925 100644
--- a/lapack-netlib/SRC/zlaswlq.f
+++ b/lapack-netlib/SRC/zlaswlq.f
@@ -1,3 +1,4 @@
+*> \brief \b ZLASWLQ
*
* Definition:
* ===========
@@ -18,9 +19,20 @@
*>
*> \verbatim
*>
-*> ZLASWLQ computes a blocked Short-Wide LQ factorization of a
-*> M-by-N matrix A, where N >= M:
-*> A = L * Q
+*> ZLASWLQ computes a blocked Tall-Skinny LQ factorization of
+*> a complexx M-by-N matrix A for M <= N:
+*>
+*> A = ( L 0 ) * Q,
+*>
+*> where:
+*>
+*> Q is a n-by-N orthogonal matrix, stored on exit in an implicit
+*> form in the elements above the digonal of the array A and in
+*> the elemenst of the array T;
+*> L is an lower-triangular M-by-M matrix stored on exit in
+*> the elements on and below the diagonal of the array A.
+*> 0 is a M-by-(N-M) zero matrix, if M < N, and is not stored.
+*>
*> \endverbatim
*
* Arguments:
@@ -150,7 +162,7 @@
SUBROUTINE ZLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
$ INFO)
*
-* -- LAPACK computational routine (version 3.7.1) --
+* -- 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. --
* June 2017
diff --git a/lapack-netlib/SRC/zlasyf_aa.f b/lapack-netlib/SRC/zlasyf_aa.f
index f321b72de..b1f1c2790 100644
--- a/lapack-netlib/SRC/zlasyf_aa.f
+++ b/lapack-netlib/SRC/zlasyf_aa.f
@@ -284,8 +284,9 @@
*
* Swap A(I1, I2+1:M) with A(I2, I2+1:M)
*
- CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
- $ A( J1+I2-1, I2+1 ), LDA )
+ IF( I2.LT.M )
+ $ CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
+ $ A( J1+I2-1, I2+1 ), LDA )
*
* Swap A(I1, I1) with A(I2,I2)
*
@@ -325,13 +326,15 @@
* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
*
- IF( A( K, J+1 ).NE.ZERO ) THEN
- ALPHA = ONE / A( K, J+1 )
- CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
- CALL ZSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
- ELSE
- CALL ZLASET( 'Full', 1, M-J-1, ZERO, ZERO,
- $ A( K, J+2 ), LDA)
+ IF( J.LT.(M-1) ) THEN
+ IF( A( K, J+1 ).NE.ZERO ) THEN
+ ALPHA = ONE / A( K, J+1 )
+ CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
+ CALL ZSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
+ ELSE
+ CALL ZLASET( 'Full', 1, M-J-1, ZERO, ZERO,
+ $ A( K, J+2 ), LDA)
+ END IF
END IF
END IF
J = J + 1
@@ -432,8 +435,9 @@
*
* Swap A(I2+1:M, I1) with A(I2+1:M, I2)
*
- CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
- $ A( I2+1, J1+I2-1 ), 1 )
+ IF( I2.LT.M )
+ $ CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
+ $ A( I2+1, J1+I2-1 ), 1 )
*
* Swap A(I1, I1) with A(I2, I2)
*
@@ -473,13 +477,15 @@
* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
*
- IF( A( J+1, K ).NE.ZERO ) THEN
- ALPHA = ONE / A( J+1, K )
- CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
- CALL ZSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
- ELSE
- CALL ZLASET( 'Full', M-J-1, 1, ZERO, ZERO,
- $ A( J+2, K ), LDA )
+ IF( J.LT.(M-1) ) THEN
+ IF( A( J+1, K ).NE.ZERO ) THEN
+ ALPHA = ONE / A( J+1, K )
+ CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
+ CALL ZSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
+ ELSE
+ CALL ZLASET( 'Full', M-J-1, 1, ZERO, ZERO,
+ $ A( J+2, K ), LDA )
+ END IF
END IF
END IF
J = J + 1
diff --git a/lapack-netlib/SRC/zlasyf_rk.f b/lapack-netlib/SRC/zlasyf_rk.f
index 664ed93f3..b6c5a27c6 100644
--- a/lapack-netlib/SRC/zlasyf_rk.f
+++ b/lapack-netlib/SRC/zlasyf_rk.f
@@ -330,7 +330,7 @@
* of A and working backwards, and compute the matrix W = U12*D
* for use in updating A11
*
-* Initilize the first entry of array E, where superdiagonal
+* Initialize the first entry of array E, where superdiagonal
* elements of D are stored
*
E( 1 ) = CZERO
@@ -658,7 +658,7 @@
* of A and working forwards, and compute the matrix W = L21*D
* for use in updating A22
*
-* Initilize the unused last entry of the subdiagonal array E.
+* Initialize the unused last entry of the subdiagonal array E.
*
E( N ) = CZERO
*
diff --git a/lapack-netlib/SRC/zlatdf.f b/lapack-netlib/SRC/zlatdf.f
index ab88570c5..4b8b5e330 100644
--- a/lapack-netlib/SRC/zlatdf.f
+++ b/lapack-netlib/SRC/zlatdf.f
@@ -261,7 +261,7 @@
*
* Solve for U- part, lockahead for RHS(N) = +-1. This is not done
* In BSOLVE and will hopefully give us a better estimate because
-* any ill-conditioning of the original matrix is transfered to U
+* any ill-conditioning of the original matrix is transferred to U
* and not to L. U(N, N) is an approximation to sigma_min(LU).
*
CALL ZCOPY( N-1, RHS, 1, WORK, 1 )
diff --git a/lapack-netlib/SRC/zlatsqr.f b/lapack-netlib/SRC/zlatsqr.f
index 1fdf3be24..0f98cae93 100644
--- a/lapack-netlib/SRC/zlatsqr.f
+++ b/lapack-netlib/SRC/zlatsqr.f
@@ -1,3 +1,4 @@
+*> \brief \b ZLATSQR
*
* Definition:
* ===========
@@ -18,9 +19,23 @@
*>
*> \verbatim
*>
-*> SLATSQR computes a blocked Tall-Skinny QR factorization of
-*> an M-by-N matrix A, where M >= N:
-*> A = Q * R .
+*> ZLATSQR computes a blocked Tall-Skinny QR factorization of
+*> a complex M-by-N matrix A for M >= N:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a M-by-M orthogonal matrix, stored on exit in an implicit
+*> form in the elements below the digonal of the array A and in
+*> the elemenst of the array T;
+*>
+*> R is an upper-triangular N-by-N matrix, stored on exit in
+*> the elements on and above the diagonal of the array A.
+*>
+*> 0 is a (M-N)-by-N zero matrix, and is not stored.
+*>
*> \endverbatim
*
* Arguments:
@@ -149,10 +164,10 @@
SUBROUTINE ZLATSQR( M, N, MB, NB, A, LDA, T, LDT, 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, MB, NB, LDT, LWORK
diff --git a/lapack-netlib/SRC/zlaunhr_col_getrfnp.f b/lapack-netlib/SRC/zlaunhr_col_getrfnp.f
new file mode 100644
index 000000000..0ab7f0349
--- /dev/null
+++ b/lapack-netlib/SRC/zlaunhr_col_getrfnp.f
@@ -0,0 +1,248 @@
+*> \brief \b ZLAUNHR_COL_GETRFNP
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZLAUNHR_COL_GETRFNP + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZLAUNHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
+*
+* .. Scalar Arguments ..
+* INTEGER INFO, LDA, M, N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 A( LDA, * ), D( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZLAUNHR_COL_GETRFNP computes the modified LU factorization without
+*> pivoting of a complex general M-by-N matrix A. The factorization has
+*> the form:
+*>
+*> A - S = L * U,
+*>
+*> where:
+*> S is a m-by-n diagonal sign matrix with the diagonal D, so that
+*> D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
+*> as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
+*> i-1 steps of Gaussian elimination. This means that the diagonal
+*> element at each step of "modified" Gaussian elimination is
+*> at least one in absolute value (so that division-by-zero not
+*> not possible during the division by the diagonal element);
+*>
+*> L is a M-by-N lower triangular matrix with unit diagonal elements
+*> (lower trapezoidal if M > N);
+*>
+*> and U is a M-by-N upper triangular matrix
+*> (upper trapezoidal if M < N).
+*>
+*> This routine is an auxiliary routine used in the Householder
+*> reconstruction routine ZUNHR_COL. In ZUNHR_COL, this routine is
+*> applied to an M-by-N matrix A with orthonormal columns, where each
+*> element is bounded by one in absolute value. With the choice of
+*> the matrix S above, one can show that the diagonal element at each
+*> step of Gaussian elimination is the largest (in absolute value) in
+*> the column on or below the diagonal, so that no pivoting is required
+*> for numerical stability [1].
+*>
+*> For more details on the Householder reconstruction algorithm,
+*> including the modified LU factorization, see [1].
+*>
+*> This is the blocked right-looking version of the algorithm,
+*> calling Level 3 BLAS to update the submatrix. To factorize a block,
+*> this routine calls the recursive routine ZLAUNHR_COL_GETRFNP2.
+*>
+*> [1] "Reconstructing Householder vectors from tall-skinny QR",
+*> G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
+*> E. Solomonik, J. Parallel Distrib. Comput.,
+*> vol. 85, pp. 3-31, 2015.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] M
+*> \verbatim
+*> M is INTEGER
+*> The number of rows of the matrix A. M >= 0.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of columns of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA,N)
+*> On entry, the M-by-N matrix to be factored.
+*> On exit, the factors L and U from the factorization
+*> A-S=L*U; the unit diagonal elements of L are not stored.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] D
+*> \verbatim
+*> D is COMPLEX*16 array, dimension min(M,N)
+*> The diagonal elements of the diagonal M-by-N sign matrix S,
+*> D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can be
+*> only ( +1.0, 0.0 ) or (-1.0, 0.0 ).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*>
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2019
+*
+*> \ingroup complex16GEcomputational
+*
+*> \par Contributors:
+* ==================
+*>
+*> \verbatim
+*>
+*> November 2019, Igor Kozachenko,
+*> Computer Science Division,
+*> University of California, Berkeley
+*>
+*> \endverbatim
+*
+* =====================================================================
+ SUBROUTINE ZLAUNHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
+ IMPLICIT NONE
+*
+* -- 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..--
+* November 2019
+*
+* .. Scalar Arguments ..
+ INTEGER INFO, LDA, M, N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 A( LDA, * ), D( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX*16 CONE
+ PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
+* ..
+* .. Local Scalars ..
+ INTEGER IINFO, J, JB, NB
+* ..
+* .. External Subroutines ..
+ EXTERNAL ZGEMM, ZLAUNHR_COL_GETRFNP2, ZTRSM, XERBLA
+* ..
+* .. External Functions ..
+ INTEGER ILAENV
+ EXTERNAL ILAENV
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX, MIN
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ IF( M.LT.0 ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+ INFO = -4
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZLAUNHR_COL_GETRFNP', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( MIN( M, N ).EQ.0 )
+ $ RETURN
+*
+* Determine the block size for this environment.
+*
+
+ NB = ILAENV( 1, 'ZLAUNHR_COL_GETRFNP', ' ', M, N, -1, -1 )
+
+ IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
+*
+* Use unblocked code.
+*
+ CALL ZLAUNHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
+ ELSE
+*
+* Use blocked code.
+*
+ DO J = 1, MIN( M, N ), NB
+ JB = MIN( MIN( M, N )-J+1, NB )
+*
+* Factor diagonal and subdiagonal blocks.
+*
+ CALL ZLAUNHR_COL_GETRFNP2( M-J+1, JB, A( J, J ), LDA,
+ $ D( J ), IINFO )
+*
+ IF( J+JB.LE.N ) THEN
+*
+* Compute block row of U.
+*
+ CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
+ $ N-J-JB+1, CONE, A( J, J ), LDA, A( J, J+JB ),
+ $ LDA )
+ IF( J+JB.LE.M ) THEN
+*
+* Update trailing submatrix.
+*
+ CALL ZGEMM( 'No transpose', 'No transpose', M-J-JB+1,
+ $ N-J-JB+1, JB, -CONE, A( J+JB, J ), LDA,
+ $ A( J, J+JB ), LDA, CONE, A( J+JB, J+JB ),
+ $ LDA )
+ END IF
+ END IF
+ END DO
+ END IF
+ RETURN
+*
+* End of ZLAUNHR_COL_GETRFNP
+*
+ END
diff --git a/lapack-netlib/SRC/zlaunhr_col_getrfnp2.f b/lapack-netlib/SRC/zlaunhr_col_getrfnp2.f
new file mode 100644
index 000000000..0057e430d
--- /dev/null
+++ b/lapack-netlib/SRC/zlaunhr_col_getrfnp2.f
@@ -0,0 +1,314 @@
+*> \brief \b ZLAUNHR_COL_GETRFNP2
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZLAUNHR_COL_GETRFNP2 + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* RECURSIVE SUBROUTINE ZLAUNHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
+*
+* .. Scalar Arguments ..
+* INTEGER INFO, LDA, M, N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 A( LDA, * ), D( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZLAUNHR_COL_GETRFNP2 computes the modified LU factorization without
+*> pivoting of a complex general M-by-N matrix A. The factorization has
+*> the form:
+*>
+*> A - S = L * U,
+*>
+*> where:
+*> S is a m-by-n diagonal sign matrix with the diagonal D, so that
+*> D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
+*> as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
+*> i-1 steps of Gaussian elimination. This means that the diagonal
+*> element at each step of "modified" Gaussian elimination is at
+*> least one in absolute value (so that division-by-zero not
+*> possible during the division by the diagonal element);
+*>
+*> L is a M-by-N lower triangular matrix with unit diagonal elements
+*> (lower trapezoidal if M > N);
+*>
+*> and U is a M-by-N upper triangular matrix
+*> (upper trapezoidal if M < N).
+*>
+*> This routine is an auxiliary routine used in the Householder
+*> reconstruction routine ZUNHR_COL. In ZUNHR_COL, this routine is
+*> applied to an M-by-N matrix A with orthonormal columns, where each
+*> element is bounded by one in absolute value. With the choice of
+*> the matrix S above, one can show that the diagonal element at each
+*> step of Gaussian elimination is the largest (in absolute value) in
+*> the column on or below the diagonal, so that no pivoting is required
+*> for numerical stability [1].
+*>
+*> For more details on the Householder reconstruction algorithm,
+*> including the modified LU factorization, see [1].
+*>
+*> This is the recursive version of the LU factorization algorithm.
+*> Denote A - S by B. The algorithm divides the matrix B into four
+*> submatrices:
+*>
+*> [ B11 | B12 ] where B11 is n1 by n1,
+*> B = [ -----|----- ] B21 is (m-n1) by n1,
+*> [ B21 | B22 ] B12 is n1 by n2,
+*> B22 is (m-n1) by n2,
+*> with n1 = min(m,n)/2, n2 = n-n1.
+*>
+*>
+*> The subroutine calls itself to factor B11, solves for B21,
+*> solves for B12, updates B22, then calls itself to factor B22.
+*>
+*> For more details on the recursive LU algorithm, see [2].
+*>
+*> ZLAUNHR_COL_GETRFNP2 is called to factorize a block by the blocked
+*> routine ZLAUNHR_COL_GETRFNP, which uses blocked code calling
+*. Level 3 BLAS to update the submatrix. However, ZLAUNHR_COL_GETRFNP2
+*> is self-sufficient and can be used without ZLAUNHR_COL_GETRFNP.
+*>
+*> [1] "Reconstructing Householder vectors from tall-skinny QR",
+*> G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
+*> E. Solomonik, J. Parallel Distrib. Comput.,
+*> vol. 85, pp. 3-31, 2015.
+*>
+*> [2] "Recursion leads to automatic variable blocking for dense linear
+*> algebra algorithms", F. Gustavson, IBM J. of Res. and Dev.,
+*> vol. 41, no. 6, pp. 737-755, 1997.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] M
+*> \verbatim
+*> M is INTEGER
+*> The number of rows of the matrix A. M >= 0.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of columns of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA,N)
+*> On entry, the M-by-N matrix to be factored.
+*> On exit, the factors L and U from the factorization
+*> A-S=L*U; the unit diagonal elements of L are not stored.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] D
+*> \verbatim
+*> D is COMPLEX*16 array, dimension min(M,N)
+*> The diagonal elements of the diagonal M-by-N sign matrix S,
+*> D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can be
+*> only ( +1.0, 0.0 ) or (-1.0, 0.0 ).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*>
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2019
+*
+*> \ingroup complex16GEcomputational
+*
+*> \par Contributors:
+* ==================
+*>
+*> \verbatim
+*>
+*> November 2019, Igor Kozachenko,
+*> Computer Science Division,
+*> University of California, Berkeley
+*>
+*> \endverbatim
+*
+* =====================================================================
+ RECURSIVE SUBROUTINE ZLAUNHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
+ IMPLICIT NONE
+*
+* -- 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..--
+* November 2019
+*
+* .. Scalar Arguments ..
+ INTEGER INFO, LDA, M, N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 A( LDA, * ), D( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ONE
+ PARAMETER ( ONE = 1.0D+0 )
+ COMPLEX*16 CONE
+ PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
+* ..
+* .. Local Scalars ..
+ DOUBLE PRECISION SFMIN
+ INTEGER I, IINFO, N1, N2
+ COMPLEX*16 Z
+* ..
+* .. External Functions ..
+ DOUBLE PRECISION DLAMCH
+ EXTERNAL DLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL ZGEMM, ZSCAL, ZTRSM, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, DBLE, DCMPLX, DIMAG, DSIGN, MAX, MIN
+* ..
+* .. Statement Functions ..
+ DOUBLE PRECISION CABS1
+* ..
+* .. Statement Function definitions ..
+ CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters
+*
+ INFO = 0
+ IF( M.LT.0 ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+ INFO = -4
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZLAUNHR_COL_GETRFNP2', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( MIN( M, N ).EQ.0 )
+ $ RETURN
+
+ IF ( M.EQ.1 ) THEN
+*
+* One row case, (also recursion termination case),
+* use unblocked code
+*
+* Transfer the sign
+*
+ D( 1 ) = DCMPLX( -DSIGN( ONE, DBLE( A( 1, 1 ) ) ) )
+*
+* Construct the row of U
+*
+ A( 1, 1 ) = A( 1, 1 ) - D( 1 )
+*
+ ELSE IF( N.EQ.1 ) THEN
+*
+* One column case, (also recursion termination case),
+* use unblocked code
+*
+* Transfer the sign
+*
+ D( 1 ) = DCMPLX( -DSIGN( ONE, DBLE( A( 1, 1 ) ) ) )
+*
+* Construct the row of U
+*
+ A( 1, 1 ) = A( 1, 1 ) - D( 1 )
+*
+* Scale the elements 2:M of the column
+*
+* Determine machine safe minimum
+*
+ SFMIN = DLAMCH('S')
+*
+* Construct the subdiagonal elements of L
+*
+ IF( CABS1( A( 1, 1 ) ) .GE. SFMIN ) THEN
+ CALL ZSCAL( M-1, CONE / A( 1, 1 ), A( 2, 1 ), 1 )
+ ELSE
+ DO I = 2, M
+ A( I, 1 ) = A( I, 1 ) / A( 1, 1 )
+ END DO
+ END IF
+*
+ ELSE
+*
+* Divide the matrix B into four submatrices
+*
+ N1 = MIN( M, N ) / 2
+ N2 = N-N1
+
+*
+* Factor B11, recursive call
+*
+ CALL ZLAUNHR_COL_GETRFNP2( N1, N1, A, LDA, D, IINFO )
+*
+* Solve for B21
+*
+ CALL ZTRSM( 'R', 'U', 'N', 'N', M-N1, N1, CONE, A, LDA,
+ $ A( N1+1, 1 ), LDA )
+*
+* Solve for B12
+*
+ CALL ZTRSM( 'L', 'L', 'N', 'U', N1, N2, CONE, A, LDA,
+ $ A( 1, N1+1 ), LDA )
+*
+* Update B22, i.e. compute the Schur complement
+* B22 := B22 - B21*B12
+*
+ CALL ZGEMM( 'N', 'N', M-N1, N2, N1, -CONE, A( N1+1, 1 ), LDA,
+ $ A( 1, N1+1 ), LDA, CONE, A( N1+1, N1+1 ), LDA )
+*
+* Factor B22, recursive call
+*
+ CALL ZLAUNHR_COL_GETRFNP2( M-N1, N2, A( N1+1, N1+1 ), LDA,
+ $ D( N1+1 ), IINFO )
+*
+ END IF
+ RETURN
+*
+* End of ZLAUNHR_COL_GETRFNP2
+*
+ END