diff --git a/lapack-netlib/SRC/ilaenv.f b/lapack-netlib/SRC/ilaenv.f
index a438ada38..7f68c383d 100644
--- a/lapack-netlib/SRC/ilaenv.f
+++ b/lapack-netlib/SRC/ilaenv.f
@@ -132,7 +132,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date November 2017
+*> \date November 2019
*
*> \ingroup OTHERauxiliary
*
@@ -162,10 +162,10 @@
* =====================================================================
INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 )
*
-* -- LAPACK auxiliary routine (version 3.8.0) --
+* -- LAPACK auxiliary 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 2017
+* November 2019
*
* .. Scalar Arguments ..
CHARACTER*( * ) NAME, OPTS
@@ -271,7 +271,16 @@
*
NB = 1
*
- IF( C2.EQ.'GE' ) THEN
+ IF( SUBNAM(2:6).EQ.'LAORH' ) THEN
+*
+* This is for *LAORHR_GETRFNP routine
+*
+ IF( SNAME ) THEN
+ NB = 32
+ ELSE
+ NB = 32
+ END IF
+ ELSE IF( C2.EQ.'GE' ) THEN
IF( C3.EQ.'TRF' ) THEN
IF( SNAME ) THEN
NB = 64
diff --git a/lapack-netlib/SRC/ilaenv2stage.f b/lapack-netlib/SRC/ilaenv2stage.f
index 3c0d34a12..db30a1b4d 100644
--- a/lapack-netlib/SRC/ilaenv2stage.f
+++ b/lapack-netlib/SRC/ilaenv2stage.f
@@ -39,9 +39,9 @@
*>
*> ILAENV2STAGE returns an INTEGER
*> if ILAENV2STAGE >= 0: ILAENV2STAGE returns the value of the parameter
-* specified by ISPEC
+*> specified by ISPEC
*> if ILAENV2STAGE < 0: if ILAENV2STAGE = -k, the k-th argument had an
-* illegal value.
+*> illegal value.
*>
*> This version provides a set of parameters which should give good,
*> but not optimal, performance on many of the currently available
diff --git a/lapack-netlib/SRC/iparam2stage.F b/lapack-netlib/SRC/iparam2stage.F
index 836e20eed..1a37300a7 100644
--- a/lapack-netlib/SRC/iparam2stage.F
+++ b/lapack-netlib/SRC/iparam2stage.F
@@ -35,7 +35,7 @@
*> \verbatim
*>
*> This program sets problem and machine dependent parameters
-*> useful for xHETRD_2STAGE, xHETRD_H@2HB, xHETRD_HB2ST,
+*> useful for xHETRD_2STAGE, xHETRD_HE2HB, xHETRD_HB2ST,
*> xGEBRD_2STAGE, xGEBRD_GE2GB, xGEBRD_GB2BD
*> and related subroutines for eigenvalue problems.
*> It is called whenever ILAENV is called with 17 <= ISPEC <= 21.
@@ -53,7 +53,7 @@
*> return.
*>
*> ISPEC=17: the optimal blocksize nb for the reduction to
-* BAND
+*> BAND
*>
*> ISPEC=18: the optimal blocksize ib for the eigenvectors
*> singular vectors update routine
@@ -90,14 +90,14 @@
*> \param[in] NBI
*> \verbatim
*> NBI is INTEGER which is the used in the reduciton,
-* (e.g., the size of the band), needed to compute workspace
-* and LHOUS2.
+*> (e.g., the size of the band), needed to compute workspace
+*> and LHOUS2.
*> \endverbatim
*>
*> \param[in] IBI
*> \verbatim
*> IBI is INTEGER which represent the IB of the reduciton,
-* needed to compute workspace and LHOUS2.
+*> needed to compute workspace and LHOUS2.
*> \endverbatim
*>
*> \param[in] NXI
diff --git a/lapack-netlib/SRC/iparmq.f b/lapack-netlib/SRC/iparmq.f
index a9212b3e0..bb711243d 100644
--- a/lapack-netlib/SRC/iparmq.f
+++ b/lapack-netlib/SRC/iparmq.f
@@ -60,7 +60,7 @@
*> invest in an (expensive) multi-shift QR sweep.
*> If the aggressive early deflation subroutine
*> finds LD converged eigenvalues from an order
-*> NW deflation window and LD.GT.(NW*NIBBLE)/100,
+*> NW deflation window and LD > (NW*NIBBLE)/100,
*> then the next QR sweep is skipped and early
*> deflation is applied immediately to the
*> remaining active diagonal block. Setting
@@ -184,8 +184,8 @@
*> This depends on ILO, IHI and NS, the
*> number of simultaneous shifts returned
*> by IPARMQ(ISPEC=15). The default for
-*> (IHI-ILO+1).LE.500 is NS. The default
-*> for (IHI-ILO+1).GT.500 is 3*NS/2.
+*> (IHI-ILO+1) <= 500 is NS. The default
+*> for (IHI-ILO+1) > 500 is 3*NS/2.
*>
*> IPARMQ(ISPEC=14) Nibble crossover point. Default: 14.
*>
diff --git a/lapack-netlib/SRC/meson.build b/lapack-netlib/SRC/meson.build
new file mode 100644
index 000000000..bad682401
--- /dev/null
+++ b/lapack-netlib/SRC/meson.build
@@ -0,0 +1,11 @@
+ALLAUX = files('ilaenv.f', 'ilaenv2stage.f', 'ieeeck.f', 'lsamen.f', 'xerbla.f', 'xerbla_array.f', 'iparmq.f', 'iparam2stage.F', 'ilaprec.f', 'ilatrans.f', 'ilauplo.f', 'iladiag.f', 'chla_transtype.f', '../INSTALL/ilaver.f', '../INSTALL/lsame.f', '../INSTALL/slamch.f')
+
+
+SCLAUX = files('sbdsdc.f', 'sbdsqr.f', 'sdisna.f', 'slabad.f', 'slacpy.f', 'sladiv.f', 'slae2.f', 'slaebz.f', 'slaed0.f', 'slaed1.f', 'slaed2.f', 'slaed3.f', 'slaed4.f', 'slaed5.f', 'slaed6.f', 'slaed7.f', 'slaed8.f', 'slaed9.f', 'slaeda.f', 'slaev2.f', 'slagtf.f', 'slagts.f', 'slamrg.f', 'slanst.f', 'slapy2.f', 'slapy3.f', 'slarnv.f', 'slarra.f', 'slarrb.f', 'slarrc.f', 'slarrd.f', 'slarre.f', 'slarrf.f', 'slarrj.f', 'slarrk.f', 'slarrr.f', 'slaneg.f', 'slartg.f', 'slaruv.f', 'slas2.f', 'slascl.f', 'slasd0.f', 'slasd1.f', 'slasd2.f', 'slasd3.f', 'slasd4.f', 'slasd5.f', 'slasd6.f', 'slasd7.f', 'slasd8.f', 'slasda.f', 'slasdq.f', 'slasdt.f', 'slaset.f', 'slasq1.f', 'slasq2.f', 'slasq3.f', 'slasq4.f', 'slasq5.f', 'slasq6.f', 'slasr.f', 'slasrt.f', 'slassq.f', 'slasv2.f', 'spttrf.f', 'sstebz.f', 'sstedc.f', 'ssteqr.f', 'ssterf.f', 'slaisnan.f', 'sisnan.f', 'slartgp.f', 'slartgs.f', '../INSTALL/second_INT_CPU_TIME.f')
+
+DZLAUX = files('dbdsdc.f', 'dbdsqr.f', 'ddisna.f', 'dlabad.f', 'dlacpy.f', 'dladiv.f', 'dlae2.f', 'dlaebz.f', 'dlaed0.f', 'dlaed1.f', 'dlaed2.f', 'dlaed3.f', 'dlaed4.f', 'dlaed5.f', 'dlaed6.f', 'dlaed7.f', 'dlaed8.f', 'dlaed9.f', 'dlaeda.f', 'dlaev2.f', 'dlagtf.f', 'dlagts.f', 'dlamrg.f', 'dlanst.f', 'dlapy2.f', 'dlapy3.f', 'dlarnv.f', 'dlarra.f', 'dlarrb.f', 'dlarrc.f', 'dlarrd.f', 'dlarre.f', 'dlarrf.f', 'dlarrj.f', 'dlarrk.f', 'dlarrr.f', 'dlaneg.f', 'dlartg.f', 'dlaruv.f', 'dlas2.f', 'dlascl.f', 'dlasd0.f', 'dlasd1.f', 'dlasd2.f', 'dlasd3.f', 'dlasd4.f', 'dlasd5.f', 'dlasd6.f', 'dlasd7.f', 'dlasd8.f', 'dlasda.f', 'dlasdq.f', 'dlasdt.f', 'dlaset.f', 'dlasq1.f', 'dlasq2.f', 'dlasq3.f', 'dlasq4.f', 'dlasq5.f', 'dlasq6.f', 'dlasr.f', 'dlasrt.f', 'dlassq.f', 'dlasv2.f', 'dpttrf.f', 'dstebz.f', 'dstedc.f', 'dsteqr.f', 'dsterf.f', 'dlaisnan.f', 'disnan.f', 'dlartgp.f', 'dlartgs.f', '../INSTALL/dlamch.f', '../INSTALL/dsecnd_INT_CPU_TIME.f')
+
+
+SLASRC = files('sbdsvdx.f', 'spotrf2.f', 'sgetrf2.f', 'sgbbrd.f', 'sgbcon.f', 'sgbequ.f', 'sgbrfs.f', 'sgbsv.f', 'sgbsvx.f', 'sgbtf2.f', 'sgbtrf.f', 'sgbtrs.f', 'sgebak.f', 'sgebal.f', 'sgebd2.f', 'sgebrd.f', 'sgecon.f', 'sgeequ.f', 'sgees.f', 'sgeesx.f', 'sgeev.f', 'sgeevx.f', 'sgehd2.f', 'sgehrd.f', 'sgelq2.f', 'sgelqf.f', 'sgels.f', 'sgelsd.f', 'sgelss.f', 'sgelsy.f', 'sgeql2.f', 'sgeqlf.f', 'sgeqp3.f', 'sgeqr2.f', 'sgeqr2p.f', 'sgeqrf.f', 'sgeqrfp.f', 'sgerfs.f', 'sgerq2.f', 'sgerqf.f', 'sgesc2.f', 'sgesdd.f', 'sgesv.f', 'sgesvd.f', 'sgesvdx.f', 'sgesvx.f', 'sgetc2.f', 'sgetf2.f', 'sgetri.f', 'sggbak.f', 'sggbal.f', 'sgges.f', 'sgges3.f', 'sggesx.f', 'sggev.f', 'sggev3.f', 'sggevx.f', 'sggglm.f', 'sgghrd.f', 'sgghd3.f', 'sgglse.f', 'sggqrf.f', 'sggrqf.f', 'sggsvd3.f', 'sggsvp3.f', 'sgtcon.f', 'sgtrfs.f', 'sgtsv.f', 'sgtsvx.f', 'sgttrf.f', 'sgttrs.f', 'sgtts2.f', 'shgeqz.f', 'shsein.f', 'shseqr.f', 'slabrd.f', 'slacon.f', 'slacn2.f', 'slaein.f', 'slaexc.f', 'slag2.f', 'slags2.f', 'slagtm.f', 'slagv2.f', 'slahqr.f', 'slahr2.f', 'slaic1.f', 'slaln2.f', 'slals0.f', 'slalsa.f', 'slalsd.f', 'slangb.f', 'slange.f', 'slangt.f', 'slanhs.f', 'slansb.f', 'slansp.f', 'slansy.f', 'slantb.f', 'slantp.f', 'slantr.f', 'slanv2.f', 'slapll.f', 'slapmt.f', 'slaqgb.f', 'slaqge.f', 'slaqp2.f', 'slaqps.f', 'slaqsb.f', 'slaqsp.f', 'slaqsy.f', 'slaqr0.f', 'slaqr1.f', 'slaqr2.f', 'slaqr3.f', 'slaqr4.f', 'slaqr5.f', 'slaqtr.f', 'slar1v.f', 'slar2v.f', 'ilaslr.f', 'ilaslc.f', 'slarf.f', 'slarfb.f', 'slarfg.f', 'slarfgp.f', 'slarft.f', 'slarfx.f', 'slarfy.f', 'slargv.f', 'slarrv.f', 'slartv.f', 'slarz.f', 'slarzb.f', 'slarzt.f', 'slaswp.f', 'slasy2.f', 'slasyf.f', 'slasyf_rook.f', 'slasyf_rk.f', 'slatbs.f', 'slatdf.f', 'slatps.f', 'slatrd.f', 'slatrs.f', 'slatrz.f', 'slauu2.f', 'slauum.f', 'sopgtr.f', 'sopmtr.f', 'sorg2l.f', 'sorg2r.f', 'sorgbr.f', 'sorghr.f', 'sorgl2.f', 'sorglq.f', 'sorgql.f', 'sorgqr.f', 'sorgr2.f', 'sorgrq.f', 'sorgtr.f', 'sorm2l.f', 'sorm2r.f', 'sorm22.f', 'sormbr.f', 'sormhr.f', 'sorml2.f', 'sormlq.f', 'sormql.f', 'sormqr.f', 'sormr2.f', 'sormr3.f', 'sormrq.f', 'sormrz.f', 'sormtr.f', 'spbcon.f', 'spbequ.f', 'spbrfs.f', 'spbstf.f', 'spbsv.f', 'spbsvx.f', 'spbtf2.f', 'spbtrf.f', 'spbtrs.f', 'spocon.f', 'spoequ.f', 'sporfs.f', 'sposv.f', 'sposvx.f', 'spotf2.f', 'spotri.f', 'spstrf.f', 'spstf2.f', 'sppcon.f', 'sppequ.f', 'spprfs.f', 'sppsv.f', 'sppsvx.f', 'spptrf.f', 'spptri.f', 'spptrs.f', 'sptcon.f', 'spteqr.f', 'sptrfs.f', 'sptsv.f', 'sptsvx.f', 'spttrs.f', 'sptts2.f', 'srscl.f', 'ssbev.f', 'ssbevd.f', 'ssbevx.f', 'ssbgst.f', 'ssbgv.f', 'ssbgvd.f', 'ssbgvx.f', 'ssbtrd.f', 'sspcon.f', 'sspev.f', 'sspevd.f', 'sspevx.f', 'sspgst.f', 'sspgv.f', 'sspgvd.f', 'sspgvx.f', 'ssprfs.f', 'sspsv.f', 'sspsvx.f', 'ssptrd.f', 'ssptrf.f', 'ssptri.f', 'ssptrs.f', 'sstegr.f', 'sstein.f', 'sstev.f', 'sstevd.f', 'sstevr.f', 'sstevx.f', 'ssycon.f', 'ssyev.f', 'ssyevd.f', 'ssyevr.f', 'ssyevx.f', 'ssygs2.f', 'ssygst.f', 'ssygv.f', 'ssygvd.f', 'ssygvx.f', 'ssyrfs.f', 'ssysv.f', 'ssysvx.f', 'ssytd2.f', 'ssytf2.f', 'ssytrd.f', 'ssytrf.f', 'ssytri.f', 'ssytri2.f', 'ssytri2x.f', 'ssyswapr.f', 'ssytrs.f', 'ssytrs2.f', 'ssyconv.f', 'ssyconvf.f', 'ssyconvf_rook.f', 'ssytf2_rook.f', 'ssytrf_rook.f', 'ssytrs_rook.f', 'ssytri_rook.f', 'ssycon_rook.f', 'ssysv_rook.f', 'ssytf2_rk.f', 'ssytrf_rk.f', 'ssytrs_3.f', 'ssytri_3.f', 'ssytri_3x.f', 'ssycon_3.f', 'ssysv_rk.f', 'slasyf_aa.f', 'ssysv_aa.f', 'ssytrf_aa.f', 'ssytrs_aa.f', 'ssysv_aa_2stage.f', 'ssytrf_aa_2stage.f', 'ssytrs_aa_2stage.f', 'stbcon.f', 'stbrfs.f', 'stbtrs.f', 'stgevc.f', 'stgex2.f', 'stgexc.f', 'stgsen.f', 'stgsja.f', 'stgsna.f', 'stgsy2.f', 'stgsyl.f', 'stpcon.f', 'stprfs.f', 'stptri.f', 'stptrs.f', 'strcon.f', 'strevc.f', 'strevc3.f', 'strexc.f', 'strrfs.f', 'strsen.f', 'strsna.f', 'strsyl.f', 'strti2.f', 'strtri.f', 'strtrs.f', 'stzrzf.f', 'sstemr.f', 'slansf.f', 'spftrf.f', 'spftri.f', 'spftrs.f', 'ssfrk.f', 'stfsm.f', 'stftri.f', 'stfttp.f', 'stfttr.f', 'stpttf.f', 'stpttr.f', 'strttf.f', 'strttp.f', 'sgejsv.f', 'sgesvj.f', 'sgsvj0.f', 'sgsvj1.f', 'sgeequb.f', 'ssyequb.f', 'spoequb.f', 'sgbequb.f', 'sbbcsd.f', 'slapmr.f', 'sorbdb.f', 'sorbdb1.f', 'sorbdb2.f', 'sorbdb3.f', 'sorbdb4.f', 'sorbdb5.f', 'sorbdb6.f', 'sorcsd.f', 'sorcsd2by1.f', 'sgeqrt.f', 'sgeqrt2.f', 'sgeqrt3.f', 'sgemqrt.f', 'stpqrt.f', 'stpqrt2.f', 'stpmqrt.f', 'stprfb.f', 'sgelqt.f', 'sgelqt3.f', 'sgemlqt.f', 'sgetsls.f', 'sgeqr.f', 'slatsqr.f', 'slamtsqr.f', 'sgemqr.f', 'sgelq.f', 'slaswlq.f', 'slamswlq.f', 'sgemlq.f', 'stplqt.f', 'stplqt2.f', 'stpmlqt.f', 'ssytrd_2stage.f', 'ssytrd_sy2sb.f', 'ssytrd_sb2st.F', 'ssb2st_kernels.f', 'ssyevd_2stage.f', 'ssyev_2stage.f', 'ssyevx_2stage.f', 'ssyevr_2stage.f', 'ssbev_2stage.f', 'ssbevx_2stage.f', 'ssbevd_2stage.f', 'ssygv_2stage.f', 'sgesvdq.f', 'scombssq.f')
+
+DSLASRC = files('spotrs.f', 'sgetrs.f', 'spotrf.f', 'sgetrf.f')
diff --git a/lapack-netlib/SRC/sbdsvdx.f b/lapack-netlib/SRC/sbdsvdx.f
index a4b1887b2..c46674c47 100644
--- a/lapack-netlib/SRC/sbdsvdx.f
+++ b/lapack-netlib/SRC/sbdsvdx.f
@@ -165,7 +165,7 @@
*>
*> \param[out] Z
*> \verbatim
-*> Z is REAL array, dimension (2*N,K) )
+*> Z is REAL array, dimension (2*N,K)
*> If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
*> contain the singular vectors of the matrix B corresponding to
*> the selected singular values, with U in rows 1 to N and V
diff --git a/lapack-netlib/SRC/scombssq.f b/lapack-netlib/SRC/scombssq.f
new file mode 100644
index 000000000..76bc0e320
--- /dev/null
+++ b/lapack-netlib/SRC/scombssq.f
@@ -0,0 +1,92 @@
+*> \brief \b SCOMBSSQ adds two scaled sum of squares quantities
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE SCOMBSSQ( V1, V2 )
+*
+* .. Array Arguments ..
+* REAL V1( 2 ), V2( 2 )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> SCOMBSSQ adds two scaled sum of squares quantities, V1 := V1 + V2.
+*> That is,
+*>
+*> V1_scale**2 * V1_sumsq := V1_scale**2 * V1_sumsq
+*> + V2_scale**2 * V2_sumsq
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in,out] V1
+*> \verbatim
+*> V1 is REAL array, dimension (2).
+*> The first scaled sum.
+*> V1(1) = V1_scale, V1(2) = V1_sumsq.
+*> \endverbatim
+*>
+*> \param[in] V2
+*> \verbatim
+*> V2 is REAL array, dimension (2).
+*> The second scaled sum.
+*> V2(1) = V2_scale, V2(2) = V2_sumsq.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2018
+*
+*> \ingroup OTHERauxiliary
+*
+* =====================================================================
+ SUBROUTINE SCOMBSSQ( V1, V2 )
+*
+* -- LAPACK auxiliary routine (version 3.7.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2018
+*
+* .. Array Arguments ..
+ REAL V1( 2 ), V2( 2 )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ZERO
+ PARAMETER ( ZERO = 0.0D+0 )
+* ..
+* .. Executable Statements ..
+*
+ IF( V1( 1 ).GE.V2( 1 ) ) THEN
+ IF( V1( 1 ).NE.ZERO ) THEN
+ V1( 2 ) = V1( 2 ) + ( V2( 1 ) / V1( 1 ) )**2 * V2( 2 )
+ END IF
+ ELSE
+ V1( 2 ) = V2( 2 ) + ( V1( 1 ) / V2( 1 ) )**2 * V1( 2 )
+ V1( 1 ) = V2( 1 )
+ END IF
+ RETURN
+*
+* End of SCOMBSSQ
+*
+ END
diff --git a/lapack-netlib/SRC/sgbrfsx.f b/lapack-netlib/SRC/sgbrfsx.f
index 032b78b80..78ae584e1 100644
--- a/lapack-netlib/SRC/sgbrfsx.f
+++ b/lapack-netlib/SRC/sgbrfsx.f
@@ -308,7 +308,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@@ -344,14 +344,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> NPARAMS is INTEGER
-*> Specifies the number of parameters set in PARAMS. If .LE. 0, the
+*> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> PARAMS is REAL array, dimension NPARAMS
-*> Specifies algorithm parameters. If an entry is .LT. 0.0, then
+*> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.
@@ -359,9 +359,9 @@
*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
*> refinement or not.
*> Default: 1.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/sgbsvxx.f b/lapack-netlib/SRC/sgbsvxx.f
index b2132325e..3c3d737b3 100644
--- a/lapack-netlib/SRC/sgbsvxx.f
+++ b/lapack-netlib/SRC/sgbsvxx.f
@@ -431,7 +431,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@@ -467,14 +467,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> NPARAMS is INTEGER
-*> Specifies the number of parameters set in PARAMS. If .LE. 0, the
+*> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> PARAMS is REAL array, dimension NPARAMS
-*> Specifies algorithm parameters. If an entry is .LT. 0.0, then
+*> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.
@@ -482,9 +482,9 @@
*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
*> refinement or not.
*> Default: 1.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/sgebak.f b/lapack-netlib/SRC/sgebak.f
index ec58bf335..5c64c8b97 100644
--- a/lapack-netlib/SRC/sgebak.f
+++ b/lapack-netlib/SRC/sgebak.f
@@ -47,10 +47,10 @@
*> \verbatim
*> JOB is CHARACTER*1
*> Specifies the type of backward transformation required:
-*> = 'N', do nothing, return immediately;
-*> = 'P', do backward transformation for permutation only;
-*> = 'S', do backward transformation for scaling only;
-*> = 'B', do backward transformations for both permutation and
+*> = 'N': do nothing, return immediately;
+*> = 'P': do backward transformation for permutation only;
+*> = 'S': do backward transformation for scaling only;
+*> = 'B': do backward transformations for both permutation and
*> scaling.
*> JOB must be the same as the argument JOB supplied to SGEBAL.
*> \endverbatim
diff --git a/lapack-netlib/SRC/sgeesx.f b/lapack-netlib/SRC/sgeesx.f
index c90de9b81..5ffa3bc37 100644
--- a/lapack-netlib/SRC/sgeesx.f
+++ b/lapack-netlib/SRC/sgeesx.f
@@ -583,7 +583,9 @@
IF( N.GT.I+1 )
$ CALL SSWAP( N-I-1, A( I, I+2 ), LDA,
$ A( I+1, I+2 ), LDA )
- CALL SSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
+ IF( WANTVS ) THEN
+ CALL SSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
+ END IF
A( I, I+1 ) = A( I+1, I )
A( I+1, I ) = ZERO
END IF
diff --git a/lapack-netlib/SRC/sgejsv.f b/lapack-netlib/SRC/sgejsv.f
index e4cbe8d0e..4ad316d99 100644
--- a/lapack-netlib/SRC/sgejsv.f
+++ b/lapack-netlib/SRC/sgejsv.f
@@ -82,7 +82,7 @@
*> desirable, then this option is advisable. The input matrix A
*> is preprocessed with QR factorization with FULL (row and
*> column) pivoting.
-*> = 'G' Computation as with 'F' with an additional estimate of the
+*> = 'G': Computation as with 'F' with an additional estimate of the
*> condition number of B, where A=D*B. If A has heavily weighted
*> rows, then using this condition number gives too pessimistic
*> error bound.
@@ -133,7 +133,7 @@
*> specified range. If A .NE. 0 is scaled so that the largest singular
*> value of c*A is around SQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
*> the licence to kill columns of A whose norm in c*A is less than
-*> SQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
+*> SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
*> where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
*> = 'N': Do not kill small columns of c*A. This option assumes that
*> BLAS and QR factorizations and triangular solvers are
@@ -230,7 +230,7 @@
*> If JOBU = 'F', then U contains on exit the M-by-M matrix of
*> the left singular vectors, including an ONB
*> of the orthogonal complement of the Range(A).
-*> If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
+*> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
*> then U is used as workspace if the procedure
*> replaces A with A^t. In that case, [V] is computed
*> in U as left singular vectors of A^t and then
@@ -252,7 +252,7 @@
*> V is REAL array, dimension ( LDV, N )
*> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
*> the right singular vectors;
-*> If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
+*> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
*> then V is used as workspace if the pprocedure
*> replaces A with A^t. In that case, [U] is computed
*> in V as right singular vectors of A^t and then
@@ -278,7 +278,7 @@
*> of A. (See the description of SVA().)
*> WORK(2) = See the description of WORK(1).
*> WORK(3) = SCONDA is an estimate for the condition number of
-*> column equilibrated A. (If JOBA .EQ. 'E' or 'G')
+*> column equilibrated A. (If JOBA = 'E' or 'G')
*> SCONDA is an estimate of SQRT(||(R^t * R)^(-1)||_1).
*> It is computed using SPOCON. It holds
*> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
@@ -297,7 +297,7 @@
*> triangular factor in the first QR factorization.
*> WORK(5) = an estimate of the scaled condition number of the
*> triangular factor in the second QR factorization.
-*> The following two parameters are computed if JOBT .EQ. 'T'.
+*> The following two parameters are computed if JOBT = 'T'.
*> They are provided for a developer/implementer who is familiar
*> with the details of the method.
*>
@@ -313,8 +313,8 @@
*> Length of WORK to confirm proper allocation of work space.
*> LWORK depends on the job:
*>
-*> If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
-*> -> .. no scaled condition estimate required (JOBE.EQ.'N'):
+*> If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and
+*> -> .. no scaled condition estimate required (JOBE = 'N'):
*> LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
*> ->> For optimal performance (blocked code) the optimal value
*> is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
@@ -330,7 +330,7 @@
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),
*> N+N*N+LWORK(DPOCON),7).
*>
-*> If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
+*> If SIGMA and the right singular vectors are needed (JOBV = 'V'),
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
*> -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
*> where NB is the optimal block size for DGEQP3, DGEQRF, DGELQ,
@@ -341,19 +341,19 @@
*> If SIGMA and the left singular vectors are needed
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
*> -> For optimal performance:
-*> if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
-*> if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
+*> if JOBU = 'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
+*> if JOBU = 'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
*> where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
*> In general, the optimal length LWORK is computed as
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
*> 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).
-*> Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or
-*> M*NB (for JOBU.EQ.'F').
+*> Here LWORK(DORMQR) equals N*NB (for JOBU = 'U') or
+*> M*NB (for JOBU = 'F').
*>
-*> If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and
-*> -> if JOBV.EQ.'V'
+*> If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
+*> -> if JOBV = 'V'
*> the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).
-*> -> if JOBV.EQ.'J' the minimal requirement is
+*> -> if JOBV = 'J' the minimal requirement is
*> LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
*> -> For optimal performance, LWORK should be additionally
*> larger than N+M*NB, where NB is the optimal block size
@@ -369,7 +369,7 @@
*> of JOBA and JOBR.
*> IWORK(2) = the number of the computed nonzero singular values
*> IWORK(3) = if nonzero, a warning message:
-*> If IWORK(3).EQ.1 then some of the column norms of A
+*> If IWORK(3) = 1 then some of the column norms of A
*> were denormalized floats. The requested high accuracy
*> is not warranted by the data.
*> \endverbatim
@@ -377,10 +377,10 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value.
-*> = 0 : successful exit;
-*> > 0 : SGEJSV did not converge in the maximal allowed number
-*> of sweeps. The computed values may be inaccurate.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value.
+*> = 0: successful exit;
+*> > 0: SGEJSV did not converge in the maximal allowed number
+*> of sweeps. The computed values may be inaccurate.
*> \endverbatim
*
* Authors:
@@ -953,7 +953,7 @@
IF ( L2ABER ) THEN
* Standard absolute error bound suffices. All sigma_i with
* sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
-* agressive enforcement of lower numerical rank by introducing a
+* aggressive enforcement of lower numerical rank by introducing a
* backward error of the order of N*EPSLN*||A||.
TEMP1 = SQRT(FLOAT(N))*EPSLN
DO 3001 p = 2, N
@@ -965,7 +965,7 @@
3001 CONTINUE
3002 CONTINUE
ELSE IF ( L2RANK ) THEN
-* .. similarly as above, only slightly more gentle (less agressive).
+* .. similarly as above, only slightly more gentle (less aggressive).
* Sudden drop on the diagonal of R1 is used as the criterion for
* close-to-rank-deficient.
TEMP1 = SQRT(SFMIN)
@@ -1294,7 +1294,7 @@
CALL SPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
$ WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
CONDR1 = ONE / SQRT(TEMP1)
-* .. here need a second oppinion on the condition number
+* .. here need a second opinion on the condition number
* .. then assume worst case scenario
* R1 is OK for inverse <=> CONDR1 .LT. FLOAT(N)
* more conservative <=> CONDR1 .LT. SQRT(FLOAT(N))
@@ -1335,7 +1335,7 @@
ELSE
*
* .. ill-conditioned case: second QRF with pivoting
-* Note that windowed pivoting would be equaly good
+* Note that windowed pivoting would be equally good
* numerically, and more run-time efficient. So, in
* an optimal implementation, the next call to SGEQP3
* should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
@@ -1388,7 +1388,7 @@
*
IF ( CONDR2 .GE. COND_OK ) THEN
* .. save the Householder vectors used for Q3
-* (this overwrittes the copy of R2, as it will not be
+* (this overwrites the copy of R2, as it will not be
* needed in this branch, but it does not overwritte the
* Huseholder vectors of Q2.).
CALL SLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
@@ -1638,7 +1638,7 @@
*
* This branch deploys a preconditioned Jacobi SVD with explicitly
* accumulated rotations. It is included as optional, mainly for
-* experimental purposes. It does perfom well, and can also be used.
+* experimental purposes. It does perform well, and can also be used.
* In this implementation, this branch will be automatically activated
* if the condition number sigma_max(A) / sigma_min(A) is predicted
* to be greater than the overflow threshold. This is because the
diff --git a/lapack-netlib/SRC/sgelq.f b/lapack-netlib/SRC/sgelq.f
index 4fe4d191d..96c4097e8 100644
--- a/lapack-netlib/SRC/sgelq.f
+++ b/lapack-netlib/SRC/sgelq.f
@@ -1,3 +1,4 @@
+*> \brief \b SGELQ
*
* Definition:
* ===========
@@ -17,7 +18,17 @@
* =============
*>
*> \verbatim
-*> SGELQ computes a LQ factorization of an M-by-N matrix A.
+*>
+*> SGELQ computes an LQ factorization of a real M-by-N matrix A:
+*>
+*> A = ( L 0 ) * Q
+*>
+*> where:
+*>
+*> Q is a N-by-N orthogonal matrix;
+*> L is an lower-triangular M-by-M matrix;
+*> 0 is a M-by-(N-M) zero matrix, if M < N.
+*>
*> \endverbatim
*
* Arguments:
@@ -138,7 +149,7 @@
*> \verbatim
*>
*> These details are particular for this LAPACK implementation. Users should not
-*> take them for granted. These details may change in the future, and are unlikely not
+*> take them for granted. These details may change in the future, and are not likely
*> true for another LAPACK implementation. These details are relevant if one wants
*> to try to understand the code. They are not part of the interface.
*>
@@ -159,10 +170,10 @@
SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N, TSIZE, LWORK
diff --git a/lapack-netlib/SRC/sgelq2.f b/lapack-netlib/SRC/sgelq2.f
index 5b1ad215a..df8128d53 100644
--- a/lapack-netlib/SRC/sgelq2.f
+++ b/lapack-netlib/SRC/sgelq2.f
@@ -33,8 +33,16 @@
*>
*> \verbatim
*>
-*> SGELQ2 computes an LQ factorization of a real m by n matrix A:
-*> A = L * Q.
+*> SGELQ2 computes an LQ factorization of a real m-by-n matrix A:
+*>
+*> A = ( L 0 ) * Q
+*>
+*> where:
+*>
+*> Q is a n-by-n orthogonal matrix;
+*> L is an lower-triangular m-by-m matrix;
+*> 0 is a m-by-(n-m) zero matrix, if m < n.
+*>
*> \endverbatim
*
* Arguments:
@@ -96,7 +104,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup realGEcomputational
*
@@ -121,10 +129,10 @@
* =====================================================================
SUBROUTINE SGELQ2( M, N, A, LDA, TAU, WORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
diff --git a/lapack-netlib/SRC/sgelqf.f b/lapack-netlib/SRC/sgelqf.f
index 99c03c0a3..90357e623 100644
--- a/lapack-netlib/SRC/sgelqf.f
+++ b/lapack-netlib/SRC/sgelqf.f
@@ -34,7 +34,15 @@
*> \verbatim
*>
*> SGELQF computes an LQ factorization of a real M-by-N matrix A:
-*> A = L * Q.
+*>
+*> A = ( L 0 ) * Q
+*>
+*> where:
+*>
+*> Q is a N-by-N orthogonal matrix;
+*> L is an lower-triangular M-by-M matrix;
+*> 0 is a M-by-(N-M) zero matrix, if M < N.
+*>
*> \endverbatim
*
* Arguments:
@@ -110,7 +118,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup realGEcomputational
*
@@ -135,10 +143,10 @@
* =====================================================================
SUBROUTINE SGELQF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, M, N
diff --git a/lapack-netlib/SRC/sgelqt.f b/lapack-netlib/SRC/sgelqt.f
index 9a93af332..64d46025c 100644
--- a/lapack-netlib/SRC/sgelqt.f
+++ b/lapack-netlib/SRC/sgelqt.f
@@ -1,3 +1,5 @@
+*> \brief \b SGELQT
+*
* Definition:
* ===========
*
diff --git a/lapack-netlib/SRC/sgelqt3.f b/lapack-netlib/SRC/sgelqt3.f
index 292ae88a3..edf5d6d30 100644
--- a/lapack-netlib/SRC/sgelqt3.f
+++ b/lapack-netlib/SRC/sgelqt3.f
@@ -1,3 +1,5 @@
+*> \brief \b SGELQT3
+*
* Definition:
* ===========
*
diff --git a/lapack-netlib/SRC/sgemlq.f b/lapack-netlib/SRC/sgemlq.f
index dedbe7752..5f2e02a8e 100644
--- a/lapack-netlib/SRC/sgemlq.f
+++ b/lapack-netlib/SRC/sgemlq.f
@@ -1,3 +1,4 @@
+*> \brief \b SGEMLQ
*
* Definition:
* ===========
@@ -143,7 +144,7 @@
*> \verbatim
*>
*> These details are particular for this LAPACK implementation. Users should not
-*> take them for granted. These details may change in the future, and are unlikely not
+*> take them for granted. These details may change in the future, and are not likely
*> true for another LAPACK implementation. These details are relevant if one wants
*> to try to understand the code. They are not part of the interface.
*>
diff --git a/lapack-netlib/SRC/sgemlqt.f b/lapack-netlib/SRC/sgemlqt.f
index a8f022bdc..37850fdf5 100644
--- a/lapack-netlib/SRC/sgemlqt.f
+++ b/lapack-netlib/SRC/sgemlqt.f
@@ -1,3 +1,5 @@
+*> \brief \b SGEMLQT
+*
* Definition:
* ===========
*
diff --git a/lapack-netlib/SRC/sgemqr.f b/lapack-netlib/SRC/sgemqr.f
index 307fc8ca9..66c5117c9 100644
--- a/lapack-netlib/SRC/sgemqr.f
+++ b/lapack-netlib/SRC/sgemqr.f
@@ -1,3 +1,4 @@
+*> \brief \b SGEMQR
*
* Definition:
* ===========
@@ -144,7 +145,7 @@
*> \verbatim
*>
*> These details are particular for this LAPACK implementation. Users should not
-*> take them for granted. These details may change in the future, and are unlikely not
+*> take them for granted. These details may change in the future, and are not likely
*> true for another LAPACK implementation. These details are relevant if one wants
*> to try to understand the code. They are not part of the interface.
*>
diff --git a/lapack-netlib/SRC/sgeqr.f b/lapack-netlib/SRC/sgeqr.f
index f939abd9d..4a6bb9ea5 100644
--- a/lapack-netlib/SRC/sgeqr.f
+++ b/lapack-netlib/SRC/sgeqr.f
@@ -1,3 +1,4 @@
+*> \brief \b SGEQR
*
* Definition:
* ===========
@@ -17,7 +18,18 @@
* =============
*>
*> \verbatim
-*> SGEQR computes a QR factorization of an M-by-N matrix A.
+*>
+*> SGEQR computes a QR factorization of a real M-by-N matrix A:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a M-by-M orthogonal matrix;
+*> R is an upper-triangular N-by-N matrix;
+*> 0 is a (M-N)-by-N zero matrix, if M > N.
+*>
*> \endverbatim
*
* Arguments:
@@ -138,7 +150,7 @@
*> \verbatim
*>
*> These details are particular for this LAPACK implementation. Users should not
-*> take them for granted. These details may change in the future, and are unlikely not
+*> take them for granted. These details may change in the future, and are not likely
*> true for another LAPACK implementation. These details are relevant if one wants
*> to try to understand the code. They are not part of the interface.
*>
@@ -160,10 +172,10 @@
SUBROUTINE SGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
$ INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N, TSIZE, LWORK
diff --git a/lapack-netlib/SRC/sgeqr2.f b/lapack-netlib/SRC/sgeqr2.f
index 3b990f825..0a1ff304f 100644
--- a/lapack-netlib/SRC/sgeqr2.f
+++ b/lapack-netlib/SRC/sgeqr2.f
@@ -33,8 +33,17 @@
*>
*> \verbatim
*>
-*> SGEQR2 computes a QR factorization of a real m by n matrix A:
-*> A = Q * R.
+*> SGEQR2 computes a QR factorization of a real m-by-n matrix A:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a m-by-m orthogonal matrix;
+*> R is an upper-triangular n-by-n matrix;
+*> 0 is a (m-n)-by-n zero matrix, if m > n.
+*>
*> \endverbatim
*
* Arguments:
@@ -96,7 +105,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup realGEcomputational
*
@@ -121,10 +130,10 @@
* =====================================================================
SUBROUTINE SGEQR2( M, N, A, LDA, TAU, WORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
diff --git a/lapack-netlib/SRC/sgeqr2p.f b/lapack-netlib/SRC/sgeqr2p.f
index f48af9d2d..08d124797 100644
--- a/lapack-netlib/SRC/sgeqr2p.f
+++ b/lapack-netlib/SRC/sgeqr2p.f
@@ -33,8 +33,18 @@
*>
*> \verbatim
*>
-*> SGEQR2P computes a QR factorization of a real m by n matrix A:
-*> A = Q * R. The diagonal entries of R are nonnegative.
+*> SGEQR2P computes a QR factorization of a real m-by-n matrix A:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a m-by-m orthogonal matrix;
+*> R is an upper-triangular n-by-n matrix with nonnegative diagonal
+*> entries;
+*> 0 is a (m-n)-by-n zero matrix, if m > n.
+*>
*> \endverbatim
*
* Arguments:
@@ -97,7 +107,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup realGEcomputational
*
@@ -124,10 +134,10 @@
* =====================================================================
SUBROUTINE SGEQR2P( M, N, A, LDA, TAU, WORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
diff --git a/lapack-netlib/SRC/sgeqrf.f b/lapack-netlib/SRC/sgeqrf.f
index 0f79c2ca5..7df495e04 100644
--- a/lapack-netlib/SRC/sgeqrf.f
+++ b/lapack-netlib/SRC/sgeqrf.f
@@ -34,7 +34,16 @@
*> \verbatim
*>
*> SGEQRF computes a QR factorization of a real M-by-N matrix A:
-*> A = Q * R.
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a M-by-M orthogonal matrix;
+*> R is an upper-triangular N-by-N matrix;
+*> 0 is a (M-N)-by-N zero matrix, if M > N.
+*>
*> \endverbatim
*
* Arguments:
@@ -111,7 +120,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup realGEcomputational
*
@@ -136,10 +145,10 @@
* =====================================================================
SUBROUTINE SGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, M, N
diff --git a/lapack-netlib/SRC/sgeqrfp.f b/lapack-netlib/SRC/sgeqrfp.f
index 654c0a13a..7f6741570 100644
--- a/lapack-netlib/SRC/sgeqrfp.f
+++ b/lapack-netlib/SRC/sgeqrfp.f
@@ -33,8 +33,18 @@
*>
*> \verbatim
*>
-*> SGEQRFP computes a QR factorization of a real M-by-N matrix A:
-*> A = Q * R. The diagonal entries of R are nonnegative.
+*> SGEQR2P computes a QR factorization of a real M-by-N matrix A:
+*>
+*> A = Q * ( R ),
+*> ( 0 )
+*>
+*> where:
+*>
+*> Q is a M-by-M orthogonal matrix;
+*> R is an upper-triangular N-by-N matrix with nonnegative diagonal
+*> entries;
+*> 0 is a (M-N)-by-N zero matrix, if M > N.
+*>
*> \endverbatim
*
* Arguments:
@@ -112,7 +122,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
-*> \date December 2016
+*> \date November 2019
*
*> \ingroup realGEcomputational
*
@@ -139,10 +149,10 @@
* =====================================================================
SUBROUTINE SGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO )
*
-* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* December 2016
+* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LWORK, M, N
diff --git a/lapack-netlib/SRC/sgerfsx.f b/lapack-netlib/SRC/sgerfsx.f
index 3f518899e..b1a1eb13d 100644
--- a/lapack-netlib/SRC/sgerfsx.f
+++ b/lapack-netlib/SRC/sgerfsx.f
@@ -283,7 +283,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@@ -319,14 +319,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> NPARAMS is INTEGER
-*> Specifies the number of parameters set in PARAMS. If .LE. 0, the
+*> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> PARAMS is REAL array, dimension NPARAMS
-*> Specifies algorithm parameters. If an entry is .LT. 0.0, then
+*> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.
@@ -334,9 +334,9 @@
*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
*> refinement or not.
*> Default: 1.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/sgesc2.f b/lapack-netlib/SRC/sgesc2.f
index c78daa334..3a6f34584 100644
--- a/lapack-netlib/SRC/sgesc2.f
+++ b/lapack-netlib/SRC/sgesc2.f
@@ -90,7 +90,7 @@
*> \verbatim
*> SCALE is REAL
*> On exit, SCALE contains the scale factor. SCALE is chosen
-*> 0 <= SCALE <= 1 to prevent owerflow in the solution.
+*> 0 <= SCALE <= 1 to prevent overflow in the solution.
*> \endverbatim
*
* Authors:
@@ -151,7 +151,7 @@
* ..
* .. Executable Statements ..
*
-* Set constant to control owerflow
+* Set constant to control overflow
*
EPS = SLAMCH( 'P' )
SMLNUM = SLAMCH( 'S' ) / EPS
diff --git a/lapack-netlib/SRC/sgesdd.f b/lapack-netlib/SRC/sgesdd.f
index 0ba2a78c7..689494dd1 100644
--- a/lapack-netlib/SRC/sgesdd.f
+++ b/lapack-netlib/SRC/sgesdd.f
@@ -322,7 +322,7 @@
*
IF( WNTQN ) THEN
* sbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
-* keep 7*N for backwards compatability.
+* keep 7*N for backwards compatibility.
BDSPAC = 7*N
ELSE
BDSPAC = 3*N*N + 4*N
@@ -448,7 +448,7 @@
*
IF( WNTQN ) THEN
* sbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6)
-* keep 7*N for backwards compatability.
+* keep 7*N for backwards compatibility.
BDSPAC = 7*M
ELSE
BDSPAC = 3*M*M + 4*M
diff --git a/lapack-netlib/SRC/sgesvdq.f b/lapack-netlib/SRC/sgesvdq.f
new file mode 100644
index 000000000..73e34862f
--- /dev/null
+++ b/lapack-netlib/SRC/sgesvdq.f
@@ -0,0 +1,1388 @@
+*> \brief SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download SGESVDQ + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE SGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
+* S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
+* WORK, LWORK, RWORK, LRWORK, INFO )
+*
+* .. Scalar Arguments ..
+* IMPLICIT NONE
+* CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
+* INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK,
+* INFO
+* ..
+* .. Array Arguments ..
+* REAL A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * )
+* REAL S( * ), RWORK( * )
+* INTEGER IWORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> SGESVDQ computes the singular value decomposition (SVD) of a real
+*> M-by-N matrix A, where M >= N. The SVD of A is written as
+*> [++] [xx] [x0] [xx]
+*> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx]
+*> [++] [xx]
+*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
+*> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
+*> of SIGMA are the singular values of A. The columns of U and V are the
+*> left and the right singular vectors of A, respectively.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] JOBA
+*> \verbatim
+*> JOBA is CHARACTER*1
+*> Specifies the level of accuracy in the computed SVD
+*> = 'A' The requested accuracy corresponds to having the backward
+*> error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F,
+*> where EPS = SLAMCH('Epsilon'). This authorises CGESVDQ to
+*> truncate the computed triangular factor in a rank revealing
+*> QR factorization whenever the truncated part is below the
+*> threshold of the order of EPS * ||A||_F. This is aggressive
+*> truncation level.
+*> = 'M' Similarly as with 'A', but the truncation is more gentle: it
+*> is allowed only when there is a drop on the diagonal of the
+*> triangular factor in the QR factorization. This is medium
+*> truncation level.
+*> = 'H' High accuracy requested. No numerical rank determination based
+*> on the rank revealing QR factorization is attempted.
+*> = 'E' Same as 'H', and in addition the condition number of column
+*> scaled A is estimated and returned in RWORK(1).
+*> N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
+*> \endverbatim
+*>
+*> \param[in] JOBP
+*> \verbatim
+*> JOBP is CHARACTER*1
+*> = 'P' The rows of A are ordered in decreasing order with respect to
+*> ||A(i,:)||_\infty. This enhances numerical accuracy at the cost
+*> of extra data movement. Recommended for numerical robustness.
+*> = 'N' No row pivoting.
+*> \endverbatim
+*>
+*> \param[in] JOBR
+*> \verbatim
+*> JOBR is CHARACTER*1
+*> = 'T' After the initial pivoted QR factorization, SGESVD is applied to
+*> the transposed R**T of the computed triangular factor R. This involves
+*> some extra data movement (matrix transpositions). Useful for
+*> experiments, research and development.
+*> = 'N' The triangular factor R is given as input to SGESVD. This may be
+*> preferred as it involves less data movement.
+*> \endverbatim
+*>
+*> \param[in] JOBU
+*> \verbatim
+*> JOBU is CHARACTER*1
+*> = 'A' All M left singular vectors are computed and returned in the
+*> matrix U. See the description of U.
+*> = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned
+*> in the matrix U. See the description of U.
+*> = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular
+*> vectors are computed and returned in the matrix U.
+*> = 'F' The N left singular vectors are returned in factored form as the
+*> product of the Q factor from the initial QR factorization and the
+*> N left singular vectors of (R**T , 0)**T. If row pivoting is used,
+*> then the necessary information on the row pivoting is stored in
+*> IWORK(N+1:N+M-1).
+*> = 'N' The left singular vectors are not computed.
+*> \endverbatim
+*>
+*> \param[in] JOBV
+*> \verbatim
+*> JOBV is CHARACTER*1
+*> = 'A', 'V' All N right singular vectors are computed and returned in
+*> the matrix V.
+*> = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular
+*> vectors are computed and returned in the matrix V. This option is
+*> allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal.
+*> = 'N' The right singular vectors are not computed.
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*> M is INTEGER
+*> The number of rows of the input matrix A. M >= 0.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of columns of the input matrix A. M >= N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is REAL array of dimensions LDA x N
+*> On entry, the input matrix A.
+*> On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains
+*> the Householder vectors as stored by SGEQP3. If JOBU = 'F', these Householder
+*> vectors together with WORK(1:N) can be used to restore the Q factors from
+*> the initial pivoted QR factorization of A. See the description of U.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER.
+*> The leading dimension of the array A. LDA >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] S
+*> \verbatim
+*> S is REAL array of dimension N.
+*> The singular values of A, ordered so that S(i) >= S(i+1).
+*> \endverbatim
+*>
+*> \param[out] U
+*> \verbatim
+*> U is REAL array, dimension
+*> LDU x M if JOBU = 'A'; see the description of LDU. In this case,
+*> on exit, U contains the M left singular vectors.
+*> LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this
+*> case, U contains the leading N or the leading NUMRANK left singular vectors.
+*> LDU x N if JOBU = 'F' ; see the description of LDU. In this case U
+*> contains N x N orthogonal matrix that can be used to form the left
+*> singular vectors.
+*> If JOBU = 'N', U is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDU
+*> \verbatim
+*> LDU is INTEGER.
+*> The leading dimension of the array U.
+*> If JOBU = 'A', 'S', 'U', 'R', LDU >= max(1,M).
+*> If JOBU = 'F', LDU >= max(1,N).
+*> Otherwise, LDU >= 1.
+*> \endverbatim
+*>
+*> \param[out] V
+*> \verbatim
+*> V is REAL array, dimension
+*> LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' .
+*> If JOBV = 'A', or 'V', V contains the N-by-N orthogonal matrix V**T;
+*> If JOBV = 'R', V contains the first NUMRANK rows of V**T (the right
+*> singular vectors, stored rowwise, of the NUMRANK largest singular values).
+*> If JOBV = 'N' and JOBA = 'E', V is used as a workspace.
+*> If JOBV = 'N', and JOBA.NE.'E', V is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*> LDV is INTEGER
+*> The leading dimension of the array V.
+*> If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N).
+*> Otherwise, LDV >= 1.
+*> \endverbatim
+*>
+*> \param[out] NUMRANK
+*> \verbatim
+*> NUMRANK is INTEGER
+*> NUMRANK is the numerical rank first determined after the rank
+*> revealing QR factorization, following the strategy specified by the
+*> value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK
+*> leading singular values and vectors are then requested in the call
+*> of SGESVD. The final value of NUMRANK might be further reduced if
+*> some singular values are computed as zeros.
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*> IWORK is INTEGER array, dimension (max(1, LIWORK)).
+*> On exit, IWORK(1:N) contains column pivoting permutation of the
+*> rank revealing QR factorization.
+*> If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence
+*> of row swaps used in row pivoting. These can be used to restore the
+*> left singular vectors in the case JOBU = 'F'.
+*>
+*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
+*> LIWORK(1) returns the minimal LIWORK.
+*> \endverbatim
+*>
+*> \param[in] LIWORK
+*> \verbatim
+*> LIWORK is INTEGER
+*> The dimension of the array IWORK.
+*> LIWORK >= N + M - 1, if JOBP = 'P' and JOBA .NE. 'E';
+*> LIWORK >= N if JOBP = 'N' and JOBA .NE. 'E';
+*> LIWORK >= N + M - 1 + N, if JOBP = 'P' and JOBA = 'E';
+*> LIWORK >= N + N if JOBP = 'N' and JOBA = 'E'.
+*
+*> If LIWORK = -1, then a workspace query is assumed; the routine
+*> only calculates and returns the optimal and minimal sizes
+*> for the WORK, IWORK, and RWORK arrays, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is REAL array, dimension (max(2, LWORK)), used as a workspace.
+*> On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains parameters
+*> needed to recover the Q factor from the QR factorization computed by
+*> SGEQP3.
+*>
+*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
+*> WORK(1) returns the optimal LWORK, and
+*> WORK(2) returns the minimal LWORK.
+*> \endverbatim
+*>
+*> \param[in,out] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of the array WORK. It is determined as follows:
+*> Let LWQP3 = 3*N+1, LWCON = 3*N, and let
+*> LWORQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U'
+*> { MAX( M, 1 ), if JOBU = 'A'
+*> LWSVD = MAX( 5*N, 1 )
+*> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ = MAX( N, 1 ),
+*> LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
+*> Then the minimal value of LWORK is:
+*> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed;
+*> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed,
+*> and a scaled condition estimate requested;
+*>
+*> = N + MAX( LWQP3, LWSVD, LWORQ ) if the singular values and the left
+*> singular vectors are requested;
+*> = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the singular values and the left
+*> singular vectors are requested, and also
+*> a scaled condition estimate requested;
+*>
+*> = N + MAX( LWQP3, LWSVD ) if the singular values and the right
+*> singular vectors are requested;
+*> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right
+*> singular vectors are requested, and also
+*> a scaled condition etimate requested;
+*>
+*> = N + MAX( LWQP3, LWSVD, LWORQ ) if the full SVD is requested with JOBV = 'R';
+*> independent of JOBR;
+*> = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is requested,
+*> JOBV = 'R' and, also a scaled condition
+*> estimate requested; independent of JOBR;
+*> = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
+*> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) ) if the
+*> full SVD is requested with JOBV = 'A' or 'V', and
+*> JOBR ='N'
+*> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
+*> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ ) )
+*> if the full SVD is requested with JOBV = 'A' or 'V', and
+*> JOBR ='N', and also a scaled condition number estimate
+*> requested.
+*> = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
+*> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) if the
+*> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
+*> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
+*> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
+*> if the full SVD is requested with JOBV = 'A' or 'V', and
+*> JOBR ='T', and also a scaled condition number estimate
+*> requested.
+*> Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK ).
+*>
+*> If LWORK = -1, then a workspace query is assumed; the routine
+*> only calculates and returns the optimal and minimal sizes
+*> for the WORK, IWORK, and RWORK arrays, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*> RWORK is REAL array, dimension (max(1, LRWORK)).
+*> On exit,
+*> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition
+*> number of column scaled A. If A = C * D where D is diagonal and C
+*> has unit columns in the Euclidean norm, then, assuming full column rank,
+*> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).
+*> Otherwise, RWORK(1) = -1.
+*> 2. RWORK(2) contains the number of singular values computed as
+*> exact zeros in SGESVD applied to the upper triangular or trapeziodal
+*> R (from the initial QR factorization). In case of early exit (no call to
+*> SGESVD, such as in the case of zero matrix) RWORK(2) = -1.
+*>
+*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
+*> RWORK(1) returns the minimal LRWORK.
+*> \endverbatim
+*>
+*> \param[in] LRWORK
+*> \verbatim
+*> LRWORK is INTEGER.
+*> The dimension of the array RWORK.
+*> If JOBP ='P', then LRWORK >= MAX(2, M).
+*> Otherwise, LRWORK >= 2
+*
+*> If LRWORK = -1, then a workspace query is assumed; the routine
+*> only calculates and returns the optimal and minimal sizes
+*> for the WORK, IWORK, and RWORK arrays, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit.
+*> < 0: if INFO = -i, the i-th argument had an illegal value.
+*> > 0: if SBDSQR did not converge, INFO specifies how many superdiagonals
+*> of an intermediate bidiagonal form B (computed in SGESVD) did not
+*> converge to zero.
+*> \endverbatim
+*
+*> \par Further Details:
+* ========================
+*>
+*> \verbatim
+*>
+*> 1. The data movement (matrix transpose) is coded using simple nested
+*> DO-loops because BLAS and LAPACK do not provide corresponding subroutines.
+*> Those DO-loops are easily identified in this source code - by the CONTINUE
+*> statements labeled with 11**. In an optimized version of this code, the
+*> nested DO loops should be replaced with calls to an optimized subroutine.
+*> 2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause
+*> column norm overflow. This is the minial precaution and it is left to the
+*> SVD routine (CGESVD) to do its own preemptive scaling if potential over-
+*> or underflows are detected. To avoid repeated scanning of the array A,
+*> an optimal implementation would do all necessary scaling before calling
+*> CGESVD and the scaling in CGESVD can be switched off.
+*> 3. Other comments related to code optimization are given in comments in the
+*> code, enlosed in [[double brackets]].
+*> \endverbatim
+*
+*> \par Bugs, examples and comments
+* ===========================
+*
+*> \verbatim
+*> Please report all bugs and send interesting examples and/or comments to
+*> drmac@math.hr. Thank you.
+*> \endverbatim
+*
+*> \par References
+* ===============
+*
+*> \verbatim
+*> [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
+*> Computing the SVD with High Accuracy. ACM Trans. Math. Softw.
+*> 44(1): 11:1-11:30 (2017)
+*>
+*> SIGMA library, xGESVDQ section updated February 2016.
+*> Developed and coded by Zlatko Drmac, Department of Mathematics
+*> University of Zagreb, Croatia, drmac@math.hr
+*> \endverbatim
+*
+*
+*> \par Contributors:
+* ==================
+*>
+*> \verbatim
+*> Developed and coded by Zlatko Drmac, Department of Mathematics
+*> University of Zagreb, Croatia, drmac@math.hr
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2018
+*
+*> \ingroup realGEsing
+*
+* =====================================================================
+ SUBROUTINE SGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
+ $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
+ $ WORK, LWORK, RWORK, LRWORK, INFO )
+* .. Scalar Arguments ..
+ IMPLICIT NONE
+ CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
+ INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK,
+ $ INFO
+* ..
+* .. Array Arguments ..
+ REAL A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * )
+ REAL S( * ), RWORK( * )
+ INTEGER IWORK( * )
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ZERO, ONE
+ PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
+* ..
+* .. Local Scalars ..
+ INTEGER IERR, IWOFF, NR, N1, OPTRATIO, p, q
+ INTEGER LWCON, LWQP3, LWRK_SGELQF, LWRK_SGESVD, LWRK_SGESVD2,
+ $ LWRK_SGEQP3, LWRK_SGEQRF, LWRK_SORMLQ, LWRK_SORMQR,
+ $ LWRK_SORMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWORQ,
+ $ LWORQ2, LWUNLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2,
+ $ IMINWRK, RMINWRK
+ LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV,
+ $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA,
+ $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR
+ REAL BIG, EPSLN, RTMP, SCONDA, SFMIN
+* ..
+* .. Local Arrays
+ REAL RDUMMY(1)
+* ..
+* .. External Subroutines (BLAS, LAPACK)
+ EXTERNAL SGELQF, SGEQP3, SGEQRF, SGESVD, SLACPY, SLAPMT,
+ $ SLASCL, SLASET, SLASWP, SSCAL, SPOCON, SORMLQ,
+ $ SORMQR, XERBLA
+* ..
+* .. External Functions (BLAS, LAPACK)
+ LOGICAL LSAME
+ INTEGER ISAMAX
+ REAL SLANGE, SNRM2, SLAMCH
+ EXTERNAL SLANGE, LSAME, ISAMAX, SNRM2, SLAMCH
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX, MIN, REAL, SQRT
+* ..
+* .. Executable Statements ..
+*
+* Test the input arguments
+*
+ WNTUS = LSAME( JOBU, 'S' ) .OR. LSAME( JOBU, 'U' )
+ WNTUR = LSAME( JOBU, 'R' )
+ WNTUA = LSAME( JOBU, 'A' )
+ WNTUF = LSAME( JOBU, 'F' )
+ LSVC0 = WNTUS .OR. WNTUR .OR. WNTUA
+ LSVEC = LSVC0 .OR. WNTUF
+ DNTWU = LSAME( JOBU, 'N' )
+*
+ WNTVR = LSAME( JOBV, 'R' )
+ WNTVA = LSAME( JOBV, 'A' ) .OR. LSAME( JOBV, 'V' )
+ RSVEC = WNTVR .OR. WNTVA
+ DNTWV = LSAME( JOBV, 'N' )
+*
+ ACCLA = LSAME( JOBA, 'A' )
+ ACCLM = LSAME( JOBA, 'M' )
+ CONDA = LSAME( JOBA, 'E' )
+ ACCLH = LSAME( JOBA, 'H' ) .OR. CONDA
+*
+ ROWPRM = LSAME( JOBP, 'P' )
+ RTRANS = LSAME( JOBR, 'T' )
+*
+ IF ( ROWPRM ) THEN
+ IF ( CONDA ) THEN
+ IMINWRK = MAX( 1, N + M - 1 + N )
+ ELSE
+ IMINWRK = MAX( 1, N + M - 1 )
+ END IF
+ RMINWRK = MAX( 2, M )
+ ELSE
+ IF ( CONDA ) THEN
+ IMINWRK = MAX( 1, N + N )
+ ELSE
+ IMINWRK = MAX( 1, N )
+ END IF
+ RMINWRK = 2
+ END IF
+ LQUERY = (LIWORK .EQ. -1 .OR. LWORK .EQ. -1 .OR. LRWORK .EQ. -1)
+ INFO = 0
+ IF ( .NOT. ( ACCLA .OR. ACCLM .OR. ACCLH ) ) THEN
+ INFO = -1
+ ELSE IF ( .NOT.( ROWPRM .OR. LSAME( JOBP, 'N' ) ) ) THEN
+ INFO = -2
+ ELSE IF ( .NOT.( RTRANS .OR. LSAME( JOBR, 'N' ) ) ) THEN
+ INFO = -3
+ ELSE IF ( .NOT.( LSVEC .OR. DNTWU ) ) THEN
+ INFO = -4
+ ELSE IF ( WNTUR .AND. WNTVA ) THEN
+ INFO = -5
+ ELSE IF ( .NOT.( RSVEC .OR. DNTWV )) THEN
+ INFO = -5
+ ELSE IF ( M.LT.0 ) THEN
+ INFO = -6
+ ELSE IF ( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
+ INFO = -7
+ ELSE IF ( LDA.LT.MAX( 1, M ) ) THEN
+ INFO = -9
+ ELSE IF ( LDU.LT.1 .OR. ( LSVC0 .AND. LDU.LT.M ) .OR.
+ $ ( WNTUF .AND. LDU.LT.N ) ) THEN
+ INFO = -12
+ ELSE IF ( LDV.LT.1 .OR. ( RSVEC .AND. LDV.LT.N ) .OR.
+ $ ( CONDA .AND. LDV.LT.N ) ) THEN
+ INFO = -14
+ ELSE IF ( LIWORK .LT. IMINWRK .AND. .NOT. LQUERY ) THEN
+ INFO = -17
+ END IF
+*
+*
+ IF ( INFO .EQ. 0 ) THEN
+* .. compute the minimal and the optimal workspace lengths
+* [[The expressions for computing the minimal and the optimal
+* values of LWORK are written with a lot of redundancy and
+* can be simplified. However, this detailed form is easier for
+* maintenance and modifications of the code.]]
+*
+* .. minimal workspace length for SGEQP3 of an M x N matrix
+ LWQP3 = 3 * N + 1
+* .. minimal workspace length for SORMQR to build left singular vectors
+ IF ( WNTUS .OR. WNTUR ) THEN
+ LWORQ = MAX( N , 1 )
+ ELSE IF ( WNTUA ) THEN
+ LWORQ = MAX( M , 1 )
+ END IF
+* .. minimal workspace length for SPOCON of an N x N matrix
+ LWCON = 3 * N
+* .. SGESVD of an N x N matrix
+ LWSVD = MAX( 5 * N, 1 )
+ IF ( LQUERY ) THEN
+ CALL SGEQP3( M, N, A, LDA, IWORK, RDUMMY, RDUMMY, -1,
+ $ IERR )
+ LWRK_SGEQP3 = INT( RDUMMY(1) )
+ IF ( WNTUS .OR. WNTUR ) THEN
+ CALL SORMQR( 'L', 'N', M, N, N, A, LDA, RDUMMY, U,
+ $ LDU, RDUMMY, -1, IERR )
+ LWRK_SORMQR = INT( RDUMMY(1) )
+ ELSE IF ( WNTUA ) THEN
+ CALL SORMQR( 'L', 'N', M, M, N, A, LDA, RDUMMY, U,
+ $ LDU, RDUMMY, -1, IERR )
+ LWRK_SORMQR = INT( RDUMMY(1) )
+ ELSE
+ LWRK_SORMQR = 0
+ END IF
+ END IF
+ MINWRK = 2
+ OPTWRK = 2
+ IF ( .NOT. (LSVEC .OR. RSVEC )) THEN
+* .. minimal and optimal sizes of the workspace if
+* only the singular values are requested
+ IF ( CONDA ) THEN
+ MINWRK = MAX( N+LWQP3, LWCON, LWSVD )
+ ELSE
+ MINWRK = MAX( N+LWQP3, LWSVD )
+ END IF
+ IF ( LQUERY ) THEN
+ CALL SGESVD( 'N', 'N', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_SGESVD = INT( RDUMMY(1) )
+ IF ( CONDA ) THEN
+ OPTWRK = MAX( N+LWRK_SGEQP3, N+LWCON, LWRK_SGESVD )
+ ELSE
+ OPTWRK = MAX( N+LWRK_SGEQP3, LWRK_SGESVD )
+ END IF
+ END IF
+ ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN
+* .. minimal and optimal sizes of the workspace if the
+* singular values and the left singular vectors are requested
+ IF ( CONDA ) THEN
+ MINWRK = N + MAX( LWQP3, LWCON, LWSVD, LWORQ )
+ ELSE
+ MINWRK = N + MAX( LWQP3, LWSVD, LWORQ )
+ END IF
+ IF ( LQUERY ) THEN
+ IF ( RTRANS ) THEN
+ CALL SGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ ELSE
+ CALL SGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ END IF
+ LWRK_SGESVD = INT( RDUMMY(1) )
+ IF ( CONDA ) THEN
+ OPTWRK = N + MAX( LWRK_SGEQP3, LWCON, LWRK_SGESVD,
+ $ LWRK_SORMQR )
+ ELSE
+ OPTWRK = N + MAX( LWRK_SGEQP3, LWRK_SGESVD,
+ $ LWRK_SORMQR )
+ END IF
+ END IF
+ ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN
+* .. minimal and optimal sizes of the workspace if the
+* singular values and the right singular vectors are requested
+ IF ( CONDA ) THEN
+ MINWRK = N + MAX( LWQP3, LWCON, LWSVD )
+ ELSE
+ MINWRK = N + MAX( LWQP3, LWSVD )
+ END IF
+ IF ( LQUERY ) THEN
+ IF ( RTRANS ) THEN
+ CALL SGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ ELSE
+ CALL SGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ END IF
+ LWRK_SGESVD = INT( RDUMMY(1) )
+ IF ( CONDA ) THEN
+ OPTWRK = N + MAX( LWRK_SGEQP3, LWCON, LWRK_SGESVD )
+ ELSE
+ OPTWRK = N + MAX( LWRK_SGEQP3, LWRK_SGESVD )
+ END IF
+ END IF
+ ELSE
+* .. minimal and optimal sizes of the workspace if the
+* full SVD is requested
+ IF ( RTRANS ) THEN
+ MINWRK = MAX( LWQP3, LWSVD, LWORQ )
+ IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON )
+ MINWRK = MINWRK + N
+ IF ( WNTVA ) THEN
+* .. minimal workspace length for N x N/2 SGEQRF
+ LWQRF = MAX( N/2, 1 )
+* .. minimal workspace lengt for N/2 x N/2 SGESVD
+ LWSVD2 = MAX( 5 * (N/2), 1 )
+ LWORQ2 = MAX( N, 1 )
+ MINWRK2 = MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2,
+ $ N/2+LWORQ2, LWORQ )
+ IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON )
+ MINWRK2 = N + MINWRK2
+ MINWRK = MAX( MINWRK, MINWRK2 )
+ END IF
+ ELSE
+ MINWRK = MAX( LWQP3, LWSVD, LWORQ )
+ IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON )
+ MINWRK = MINWRK + N
+ IF ( WNTVA ) THEN
+* .. minimal workspace length for N/2 x N SGELQF
+ LWLQF = MAX( N/2, 1 )
+ LWSVD2 = MAX( 5 * (N/2), 1 )
+ LWUNLQ = MAX( N , 1 )
+ MINWRK2 = MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2,
+ $ N/2+LWUNLQ, LWORQ )
+ IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON )
+ MINWRK2 = N + MINWRK2
+ MINWRK = MAX( MINWRK, MINWRK2 )
+ END IF
+ END IF
+ IF ( LQUERY ) THEN
+ IF ( RTRANS ) THEN
+ CALL SGESVD( 'O', 'A', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_SGESVD = INT( RDUMMY(1) )
+ OPTWRK = MAX(LWRK_SGEQP3,LWRK_SGESVD,LWRK_SORMQR)
+ IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON )
+ OPTWRK = N + OPTWRK
+ IF ( WNTVA ) THEN
+ CALL SGEQRF(N,N/2,U,LDU,RDUMMY,RDUMMY,-1,IERR)
+ LWRK_SGEQRF = INT( RDUMMY(1) )
+ CALL SGESVD( 'S', 'O', N/2,N/2, V,LDV, S, U,LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_SGESVD2 = INT( RDUMMY(1) )
+ CALL SORMQR( 'R', 'C', N, N, N/2, U, LDU, RDUMMY,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_SORMQR2 = INT( RDUMMY(1) )
+ OPTWRK2 = MAX( LWRK_SGEQP3, N/2+LWRK_SGEQRF,
+ $ N/2+LWRK_SGESVD2, N/2+LWRK_SORMQR2 )
+ IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON )
+ OPTWRK2 = N + OPTWRK2
+ OPTWRK = MAX( OPTWRK, OPTWRK2 )
+ END IF
+ ELSE
+ CALL SGESVD( 'S', 'O', N, N, A, LDA, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_SGESVD = INT( RDUMMY(1) )
+ OPTWRK = MAX(LWRK_SGEQP3,LWRK_SGESVD,LWRK_SORMQR)
+ IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON )
+ OPTWRK = N + OPTWRK
+ IF ( WNTVA ) THEN
+ CALL SGELQF(N/2,N,U,LDU,RDUMMY,RDUMMY,-1,IERR)
+ LWRK_SGELQF = INT( RDUMMY(1) )
+ CALL SGESVD( 'S','O', N/2,N/2, V, LDV, S, U, LDU,
+ $ V, LDV, RDUMMY, -1, IERR )
+ LWRK_SGESVD2 = INT( RDUMMY(1) )
+ CALL SORMLQ( 'R', 'N', N, N, N/2, U, LDU, RDUMMY,
+ $ V, LDV, RDUMMY,-1,IERR )
+ LWRK_SORMLQ = INT( RDUMMY(1) )
+ OPTWRK2 = MAX( LWRK_SGEQP3, N/2+LWRK_SGELQF,
+ $ N/2+LWRK_SGESVD2, N/2+LWRK_SORMLQ )
+ IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON )
+ OPTWRK2 = N + OPTWRK2
+ OPTWRK = MAX( OPTWRK, OPTWRK2 )
+ END IF
+ END IF
+ END IF
+ END IF
+*
+ MINWRK = MAX( 2, MINWRK )
+ OPTWRK = MAX( 2, OPTWRK )
+ IF ( LWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = -19
+*
+ END IF
+*
+ IF (INFO .EQ. 0 .AND. LRWORK .LT. RMINWRK .AND. .NOT. LQUERY) THEN
+ INFO = -21
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'SGESVDQ', -INFO )
+ RETURN
+ ELSE IF ( LQUERY ) THEN
+*
+* Return optimal workspace
+*
+ IWORK(1) = IMINWRK
+ WORK(1) = OPTWRK
+ WORK(2) = MINWRK
+ RWORK(1) = RMINWRK
+ RETURN
+ END IF
+*
+* Quick return if the matrix is void.
+*
+ IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) THEN
+* .. all output is void.
+ RETURN
+ END IF
+*
+ BIG = SLAMCH('O')
+ ASCALED = .FALSE.
+ IWOFF = 1
+ IF ( ROWPRM ) THEN
+ IWOFF = M
+* .. reordering the rows in decreasing sequence in the
+* ell-infinity norm - this enhances numerical robustness in
+* the case of differently scaled rows.
+ DO 1904 p = 1, M
+* RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) )
+* [[SLANGE will return NaN if an entry of the p-th row is Nan]]
+ RWORK(p) = SLANGE( 'M', 1, N, A(p,1), LDA, RDUMMY )
+* .. check for NaN's and Inf's
+ IF ( ( RWORK(p) .NE. RWORK(p) ) .OR.
+ $ ( (RWORK(p)*ZERO) .NE. ZERO ) ) THEN
+ INFO = -8
+ CALL XERBLA( 'SGESVDQ', -INFO )
+ RETURN
+ END IF
+ 1904 CONTINUE
+ DO 1952 p = 1, M - 1
+ q = ISAMAX( M-p+1, RWORK(p), 1 ) + p - 1
+ IWORK(N+p) = q
+ IF ( p .NE. q ) THEN
+ RTMP = RWORK(p)
+ RWORK(p) = RWORK(q)
+ RWORK(q) = RTMP
+ END IF
+ 1952 CONTINUE
+*
+ IF ( RWORK(1) .EQ. ZERO ) THEN
+* Quick return: A is the M x N zero matrix.
+ NUMRANK = 0
+ CALL SLASET( 'G', N, 1, ZERO, ZERO, S, N )
+ IF ( WNTUS ) CALL SLASET('G', M, N, ZERO, ONE, U, LDU)
+ IF ( WNTUA ) CALL SLASET('G', M, M, ZERO, ONE, U, LDU)
+ IF ( WNTVA ) CALL SLASET('G', N, N, ZERO, ONE, V, LDV)
+ IF ( WNTUF ) THEN
+ CALL SLASET( 'G', N, 1, ZERO, ZERO, WORK, N )
+ CALL SLASET( 'G', M, N, ZERO, ONE, U, LDU )
+ END IF
+ DO 5001 p = 1, N
+ IWORK(p) = p
+ 5001 CONTINUE
+ IF ( ROWPRM ) THEN
+ DO 5002 p = N + 1, N + M - 1
+ IWORK(p) = p - N
+ 5002 CONTINUE
+ END IF
+ IF ( CONDA ) RWORK(1) = -1
+ RWORK(2) = -1
+ RETURN
+ END IF
+*
+ IF ( RWORK(1) .GT. BIG / SQRT(REAL(M)) ) THEN
+* .. to prevent overflow in the QR factorization, scale the
+* matrix by 1/sqrt(M) if too large entry detected
+ CALL SLASCL('G',0,0,SQRT(REAL(M)),ONE, M,N, A,LDA, IERR)
+ ASCALED = .TRUE.
+ END IF
+ CALL SLASWP( N, A, LDA, 1, M-1, IWORK(N+1), 1 )
+ END IF
+*
+* .. At this stage, preemptive scaling is done only to avoid column
+* norms overflows during the QR factorization. The SVD procedure should
+* have its own scaling to save the singular values from overflows and
+* underflows. That depends on the SVD procedure.
+*
+ IF ( .NOT.ROWPRM ) THEN
+ RTMP = SLANGE( 'M', M, N, A, LDA, RDUMMY )
+ IF ( ( RTMP .NE. RTMP ) .OR.
+ $ ( (RTMP*ZERO) .NE. ZERO ) ) THEN
+ INFO = -8
+ CALL XERBLA( 'SGESVDQ', -INFO )
+ RETURN
+ END IF
+ IF ( RTMP .GT. BIG / SQRT(REAL(M)) ) THEN
+* .. to prevent overflow in the QR factorization, scale the
+* matrix by 1/sqrt(M) if too large entry detected
+ CALL SLASCL('G',0,0, SQRT(REAL(M)),ONE, M,N, A,LDA, IERR)
+ ASCALED = .TRUE.
+ END IF
+ END IF
+*
+* .. QR factorization with column pivoting
+*
+* A * P = Q * [ R ]
+* [ 0 ]
+*
+ DO 1963 p = 1, N
+* .. all columns are free columns
+ IWORK(p) = 0
+ 1963 CONTINUE
+ CALL SGEQP3( M, N, A, LDA, IWORK, WORK, WORK(N+1), LWORK-N,
+ $ IERR )
+*
+* If the user requested accuracy level allows truncation in the
+* computed upper triangular factor, the matrix R is examined and,
+* if possible, replaced with its leading upper trapezoidal part.
+*
+ EPSLN = SLAMCH('E')
+ SFMIN = SLAMCH('S')
+* SMALL = SFMIN / EPSLN
+ NR = N
+*
+ IF ( ACCLA ) THEN
+*
+* Standard absolute error bound suffices. All sigma_i with
+* sigma_i < N*EPS*||A||_F are flushed to zero. This is an
+* aggressive enforcement of lower numerical rank by introducing a
+* backward error of the order of N*EPS*||A||_F.
+ NR = 1
+ RTMP = SQRT(REAL(N))*EPSLN
+ DO 3001 p = 2, N
+ IF ( ABS(A(p,p)) .LT. (RTMP*ABS(A(1,1))) ) GO TO 3002
+ NR = NR + 1
+ 3001 CONTINUE
+ 3002 CONTINUE
+*
+ ELSEIF ( ACCLM ) THEN
+* .. similarly as above, only slightly more gentle (less aggressive).
+* Sudden drop on the diagonal of R is used as the criterion for being
+* close-to-rank-deficient. The threshold is set to EPSLN=SLAMCH('E').
+* [[This can be made more flexible by replacing this hard-coded value
+* with a user specified threshold.]] Also, the values that underflow
+* will be truncated.
+ NR = 1
+ DO 3401 p = 2, N
+ IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
+ $ ( ABS(A(p,p)) .LT. SFMIN ) ) GO TO 3402
+ NR = NR + 1
+ 3401 CONTINUE
+ 3402 CONTINUE
+*
+ ELSE
+* .. RRQR not authorized to determine numerical rank except in the
+* obvious case of zero pivots.
+* .. inspect R for exact zeros on the diagonal;
+* R(i,i)=0 => R(i:N,i:N)=0.
+ NR = 1
+ DO 3501 p = 2, N
+ IF ( ABS(A(p,p)) .EQ. ZERO ) GO TO 3502
+ NR = NR + 1
+ 3501 CONTINUE
+ 3502 CONTINUE
+*
+ IF ( CONDA ) THEN
+* Estimate the scaled condition number of A. Use the fact that it is
+* the same as the scaled condition number of R.
+* .. V is used as workspace
+ CALL SLACPY( 'U', N, N, A, LDA, V, LDV )
+* Only the leading NR x NR submatrix of the triangular factor
+* is considered. Only if NR=N will this give a reliable error
+* bound. However, even for NR < N, this can be used on an
+* expert level and obtain useful information in the sense of
+* perturbation theory.
+ DO 3053 p = 1, NR
+ RTMP = SNRM2( p, V(1,p), 1 )
+ CALL SSCAL( p, ONE/RTMP, V(1,p), 1 )
+ 3053 CONTINUE
+ IF ( .NOT. ( LSVEC .OR. RSVEC ) ) THEN
+ CALL SPOCON( 'U', NR, V, LDV, ONE, RTMP,
+ $ WORK, IWORK(N+IWOFF), IERR )
+ ELSE
+ CALL SPOCON( 'U', NR, V, LDV, ONE, RTMP,
+ $ WORK(N+1), IWORK(N+IWOFF), IERR )
+ END IF
+ SCONDA = ONE / SQRT(RTMP)
+* For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
+* N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
+* See the reference [1] for more details.
+ END IF
+*
+ ENDIF
+*
+ IF ( WNTUR ) THEN
+ N1 = NR
+ ELSE IF ( WNTUS .OR. WNTUF) THEN
+ N1 = N
+ ELSE IF ( WNTUA ) THEN
+ N1 = M
+ END IF
+*
+ IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
+*.......................................................................
+* .. only the singular values are requested
+*.......................................................................
+ IF ( RTRANS ) THEN
+*
+* .. compute the singular values of R**T = [A](1:NR,1:N)**T
+* .. set the lower triangle of [A] to [A](1:NR,1:N)**T and
+* the upper triangle of [A] to zero.
+ DO 1146 p = 1, MIN( N, NR )
+ DO 1147 q = p + 1, N
+ A(q,p) = A(p,q)
+ IF ( q .LE. NR ) A(p,q) = ZERO
+ 1147 CONTINUE
+ 1146 CONTINUE
+*
+ CALL SGESVD( 'N', 'N', N, NR, A, LDA, S, U, LDU,
+ $ V, LDV, WORK, LWORK, INFO )
+*
+ ELSE
+*
+* .. compute the singular values of R = [A](1:NR,1:N)
+*
+ IF ( NR .GT. 1 )
+ $ CALL SLASET( 'L', NR-1,NR-1, ZERO,ZERO, A(2,1), LDA )
+ CALL SGESVD( 'N', 'N', NR, N, A, LDA, S, U, LDU,
+ $ V, LDV, WORK, LWORK, INFO )
+*
+ END IF
+*
+ ELSE IF ( LSVEC .AND. ( .NOT. RSVEC) ) THEN
+*.......................................................................
+* .. the singular values and the left singular vectors requested
+*.......................................................................""""""""
+ IF ( RTRANS ) THEN
+* .. apply SGESVD to R**T
+* .. copy R**T into [U] and overwrite [U] with the right singular
+* vectors of R
+ DO 1192 p = 1, NR
+ DO 1193 q = p, N
+ U(q,p) = A(p,q)
+ 1193 CONTINUE
+ 1192 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL SLASET( 'U', NR-1,NR-1, ZERO,ZERO, U(1,2), LDU )
+* .. the left singular vectors not computed, the NR right singular
+* vectors overwrite [U](1:NR,1:NR) as transposed. These
+* will be pre-multiplied by Q to build the left singular vectors of A.
+ CALL SGESVD( 'N', 'O', N, NR, U, LDU, S, U, LDU,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+*
+ DO 1119 p = 1, NR
+ DO 1120 q = p + 1, NR
+ RTMP = U(q,p)
+ U(q,p) = U(p,q)
+ U(p,q) = RTMP
+ 1120 CONTINUE
+ 1119 CONTINUE
+*
+ ELSE
+* .. apply SGESVD to R
+* .. copy R into [U] and overwrite [U] with the left singular vectors
+ CALL SLACPY( 'U', NR, N, A, LDA, U, LDU )
+ IF ( NR .GT. 1 )
+ $ CALL SLASET( 'L', NR-1, NR-1, ZERO, ZERO, U(2,1), LDU )
+* .. the right singular vectors not computed, the NR left singular
+* vectors overwrite [U](1:NR,1:NR)
+ CALL SGESVD( 'O', 'N', NR, N, U, LDU, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+* .. now [U](1:NR,1:NR) contains the NR left singular vectors of
+* R. These will be pre-multiplied by Q to build the left singular
+* vectors of A.
+ END IF
+*
+* .. assemble the left singular vector matrix U of dimensions
+* (M x NR) or (M x N) or (M x M).
+ IF ( ( NR .LT. M ) .AND. ( .NOT.WNTUF ) ) THEN
+ CALL SLASET('A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL SLASET( 'A',NR,N1-NR,ZERO,ZERO,U(1,NR+1), LDU )
+ CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1), LDU )
+ END IF
+ END IF
+*
+* The Q matrix from the first QRF is built into the left singular
+* vectors matrix U.
+*
+ IF ( .NOT.WNTUF )
+ $ CALL SORMQR( 'L', 'N', M, N1, N, A, LDA, WORK, U,
+ $ LDU, WORK(N+1), LWORK-N, IERR )
+ IF ( ROWPRM .AND. .NOT.WNTUF )
+ $ CALL SLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 )
+*
+ ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
+*.......................................................................
+* .. the singular values and the right singular vectors requested
+*.......................................................................
+ IF ( RTRANS ) THEN
+* .. apply SGESVD to R**T
+* .. copy R**T into V and overwrite V with the left singular vectors
+ DO 1165 p = 1, NR
+ DO 1166 q = p, N
+ V(q,p) = (A(p,q))
+ 1166 CONTINUE
+ 1165 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL SLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
+* .. the left singular vectors of R**T overwrite V, the right singular
+* vectors not computed
+ IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
+ CALL SGESVD( 'O', 'N', N, NR, V, LDV, S, U, LDU,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+*
+ DO 1121 p = 1, NR
+ DO 1122 q = p + 1, NR
+ RTMP = V(q,p)
+ V(q,p) = V(p,q)
+ V(p,q) = RTMP
+ 1122 CONTINUE
+ 1121 CONTINUE
+*
+ IF ( NR .LT. N ) THEN
+ DO 1103 p = 1, NR
+ DO 1104 q = NR + 1, N
+ V(p,q) = V(q,p)
+ 1104 CONTINUE
+ 1103 CONTINUE
+ END IF
+ CALL SLAPMT( .FALSE., NR, N, V, LDV, IWORK )
+ ELSE
+* .. need all N right singular vectors and NR < N
+* [!] This is simple implementation that augments [V](1:N,1:NR)
+* by padding a zero block. In the case NR << N, a more efficient
+* way is to first use the QR factorization. For more details
+* how to implement this, see the " FULL SVD " branch.
+ CALL SLASET('G', N, N-NR, ZERO, ZERO, V(1,NR+1), LDV)
+ CALL SGESVD( 'O', 'N', N, N, V, LDV, S, U, LDU,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+*
+ DO 1123 p = 1, N
+ DO 1124 q = p + 1, N
+ RTMP = V(q,p)
+ V(q,p) = V(p,q)
+ V(p,q) = RTMP
+ 1124 CONTINUE
+ 1123 CONTINUE
+ CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK )
+ END IF
+*
+ ELSE
+* .. aply SGESVD to R
+* .. copy R into V and overwrite V with the right singular vectors
+ CALL SLACPY( 'U', NR, N, A, LDA, V, LDV )
+ IF ( NR .GT. 1 )
+ $ CALL SLASET( 'L', NR-1, NR-1, ZERO, ZERO, V(2,1), LDV )
+* .. the right singular vectors overwrite V, the NR left singular
+* vectors stored in U(1:NR,1:NR)
+ IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
+ CALL SGESVD( 'N', 'O', NR, N, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+ CALL SLAPMT( .FALSE., NR, N, V, LDV, IWORK )
+* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
+ ELSE
+* .. need all N right singular vectors and NR < N
+* [!] This is simple implementation that augments [V](1:NR,1:N)
+* by padding a zero block. In the case NR << N, a more efficient
+* way is to first use the LQ factorization. For more details
+* how to implement this, see the " FULL SVD " branch.
+ CALL SLASET('G', N-NR, N, ZERO,ZERO, V(NR+1,1), LDV)
+ CALL SGESVD( 'N', 'O', N, N, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+ CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK )
+ END IF
+* .. now [V] contains the transposed matrix of the right singular
+* vectors of A.
+ END IF
+*
+ ELSE
+*.......................................................................
+* .. FULL SVD requested
+*.......................................................................
+ IF ( RTRANS ) THEN
+*
+* .. apply SGESVD to R**T [[this option is left for R&D&T]]
+*
+ IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
+* .. copy R**T into [V] and overwrite [V] with the left singular
+* vectors of R**T
+ DO 1168 p = 1, NR
+ DO 1169 q = p, N
+ V(q,p) = A(p,q)
+ 1169 CONTINUE
+ 1168 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL SLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
+*
+* .. the left singular vectors of R**T overwrite [V], the NR right
+* singular vectors of R**T stored in [U](1:NR,1:NR) as transposed
+ CALL SGESVD( 'O', 'A', N, NR, V, LDV, S, V, LDV,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+* .. assemble V
+ DO 1115 p = 1, NR
+ DO 1116 q = p + 1, NR
+ RTMP = V(q,p)
+ V(q,p) = V(p,q)
+ V(p,q) = RTMP
+ 1116 CONTINUE
+ 1115 CONTINUE
+ IF ( NR .LT. N ) THEN
+ DO 1101 p = 1, NR
+ DO 1102 q = NR+1, N
+ V(p,q) = V(q,p)
+ 1102 CONTINUE
+ 1101 CONTINUE
+ END IF
+ CALL SLAPMT( .FALSE., NR, N, V, LDV, IWORK )
+*
+ DO 1117 p = 1, NR
+ DO 1118 q = p + 1, NR
+ RTMP = U(q,p)
+ U(q,p) = U(p,q)
+ U(p,q) = RTMP
+ 1118 CONTINUE
+ 1117 CONTINUE
+*
+ IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL SLASET('A', M-NR,NR, ZERO,ZERO, U(NR+1,1), LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
+ CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1), LDU )
+ END IF
+ END IF
+*
+ ELSE
+* .. need all N right singular vectors and NR < N
+* .. copy R**T into [V] and overwrite [V] with the left singular
+* vectors of R**T
+* [[The optimal ratio N/NR for using QRF instead of padding
+* with zeros. Here hard coded to 2; it must be at least
+* two due to work space constraints.]]
+* OPTRATIO = ILAENV(6, 'SGESVD', 'S' // 'O', NR,N,0,0)
+* OPTRATIO = MAX( OPTRATIO, 2 )
+ OPTRATIO = 2
+ IF ( OPTRATIO*NR .GT. N ) THEN
+ DO 1198 p = 1, NR
+ DO 1199 q = p, N
+ V(q,p) = A(p,q)
+ 1199 CONTINUE
+ 1198 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL SLASET('U',NR-1,NR-1, ZERO,ZERO, V(1,2),LDV)
+*
+ CALL SLASET('A',N,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
+ CALL SGESVD( 'O', 'A', N, N, V, LDV, S, V, LDV,
+ $ U, LDU, WORK(N+1), LWORK-N, INFO )
+*
+ DO 1113 p = 1, N
+ DO 1114 q = p + 1, N
+ RTMP = V(q,p)
+ V(q,p) = V(p,q)
+ V(p,q) = RTMP
+ 1114 CONTINUE
+ 1113 CONTINUE
+ CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK )
+* .. assemble the left singular vector matrix U of dimensions
+* (M x N1), i.e. (M x N) or (M x M).
+*
+ DO 1111 p = 1, N
+ DO 1112 q = p + 1, N
+ RTMP = U(q,p)
+ U(q,p) = U(p,q)
+ U(p,q) = RTMP
+ 1112 CONTINUE
+ 1111 CONTINUE
+*
+ IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL SLASET('A',M-N,N,ZERO,ZERO,U(N+1,1),LDU)
+ IF ( N .LT. N1 ) THEN
+ CALL SLASET('A',N,N1-N,ZERO,ZERO,U(1,N+1),LDU)
+ CALL SLASET('A',M-N,N1-N,ZERO,ONE,
+ $ U(N+1,N+1), LDU )
+ END IF
+ END IF
+ ELSE
+* .. copy R**T into [U] and overwrite [U] with the right
+* singular vectors of R
+ DO 1196 p = 1, NR
+ DO 1197 q = p, N
+ U(q,NR+p) = A(p,q)
+ 1197 CONTINUE
+ 1196 CONTINUE
+ IF ( NR .GT. 1 )
+ $ CALL SLASET('U',NR-1,NR-1,ZERO,ZERO,U(1,NR+2),LDU)
+ CALL SGEQRF( N, NR, U(1,NR+1), LDU, WORK(N+1),
+ $ WORK(N+NR+1), LWORK-N-NR, IERR )
+ DO 1143 p = 1, NR
+ DO 1144 q = 1, N
+ V(q,p) = U(p,NR+q)
+ 1144 CONTINUE
+ 1143 CONTINUE
+ CALL SLASET('U',NR-1,NR-1,ZERO,ZERO,V(1,2),LDV)
+ CALL SGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU,
+ $ V,LDV, WORK(N+NR+1),LWORK-N-NR, INFO )
+ CALL SLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
+ CALL SLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
+ CALL SLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
+ CALL SORMQR('R','C', N, N, NR, U(1,NR+1), LDU,
+ $ WORK(N+1),V,LDV,WORK(N+NR+1),LWORK-N-NR,IERR)
+ CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK )
+* .. assemble the left singular vector matrix U of dimensions
+* (M x NR) or (M x N) or (M x M).
+ IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL SLASET('A',M-NR,NR,ZERO,ZERO,U(NR+1,1),LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
+ CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1),LDU)
+ END IF
+ END IF
+ END IF
+ END IF
+*
+ ELSE
+*
+* .. apply SGESVD to R [[this is the recommended option]]
+*
+ IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
+* .. copy R into [V] and overwrite V with the right singular vectors
+ CALL SLACPY( 'U', NR, N, A, LDA, V, LDV )
+ IF ( NR .GT. 1 )
+ $ CALL SLASET( 'L', NR-1,NR-1, ZERO,ZERO, V(2,1), LDV )
+* .. the right singular vectors of R overwrite [V], the NR left
+* singular vectors of R stored in [U](1:NR,1:NR)
+ CALL SGESVD( 'S', 'O', NR, N, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+ CALL SLAPMT( .FALSE., NR, N, V, LDV, IWORK )
+* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
+* .. assemble the left singular vector matrix U of dimensions
+* (M x NR) or (M x N) or (M x M).
+ IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL SLASET('A', M-NR,NR, ZERO,ZERO, U(NR+1,1), LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
+ CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1), LDU )
+ END IF
+ END IF
+*
+ ELSE
+* .. need all N right singular vectors and NR < N
+* .. the requested number of the left singular vectors
+* is then N1 (N or M)
+* [[The optimal ratio N/NR for using LQ instead of padding
+* with zeros. Here hard coded to 2; it must be at least
+* two due to work space constraints.]]
+* OPTRATIO = ILAENV(6, 'SGESVD', 'S' // 'O', NR,N,0,0)
+* OPTRATIO = MAX( OPTRATIO, 2 )
+ OPTRATIO = 2
+ IF ( OPTRATIO * NR .GT. N ) THEN
+ CALL SLACPY( 'U', NR, N, A, LDA, V, LDV )
+ IF ( NR .GT. 1 )
+ $ CALL SLASET('L', NR-1,NR-1, ZERO,ZERO, V(2,1),LDV)
+* .. the right singular vectors of R overwrite [V], the NR left
+* singular vectors of R stored in [U](1:NR,1:NR)
+ CALL SLASET('A', N-NR,N, ZERO,ZERO, V(NR+1,1),LDV)
+ CALL SGESVD( 'S', 'O', N, N, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+1), LWORK-N, INFO )
+ CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK )
+* .. now [V] contains the transposed matrix of the right
+* singular vectors of A. The leading N left singular vectors
+* are in [U](1:N,1:N)
+* .. assemble the left singular vector matrix U of dimensions
+* (M x N1), i.e. (M x N) or (M x M).
+ IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL SLASET('A',M-N,N,ZERO,ZERO,U(N+1,1),LDU)
+ IF ( N .LT. N1 ) THEN
+ CALL SLASET('A',N,N1-N,ZERO,ZERO,U(1,N+1),LDU)
+ CALL SLASET( 'A',M-N,N1-N,ZERO,ONE,
+ $ U(N+1,N+1), LDU )
+ END IF
+ END IF
+ ELSE
+ CALL SLACPY( 'U', NR, N, A, LDA, U(NR+1,1), LDU )
+ IF ( NR .GT. 1 )
+ $ CALL SLASET('L',NR-1,NR-1,ZERO,ZERO,U(NR+2,1),LDU)
+ CALL SGELQF( NR, N, U(NR+1,1), LDU, WORK(N+1),
+ $ WORK(N+NR+1), LWORK-N-NR, IERR )
+ CALL SLACPY('L',NR,NR,U(NR+1,1),LDU,V,LDV)
+ IF ( NR .GT. 1 )
+ $ CALL SLASET('U',NR-1,NR-1,ZERO,ZERO,V(1,2),LDV)
+ CALL SGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU,
+ $ V, LDV, WORK(N+NR+1), LWORK-N-NR, INFO )
+ CALL SLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
+ CALL SLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
+ CALL SLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
+ CALL SORMLQ('R','N',N,N,NR,U(NR+1,1),LDU,WORK(N+1),
+ $ V, LDV, WORK(N+NR+1),LWORK-N-NR,IERR)
+ CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK )
+* .. assemble the left singular vector matrix U of dimensions
+* (M x NR) or (M x N) or (M x M).
+ IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
+ CALL SLASET('A',M-NR,NR,ZERO,ZERO,U(NR+1,1),LDU)
+ IF ( NR .LT. N1 ) THEN
+ CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
+ CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE,
+ $ U(NR+1,NR+1), LDU )
+ END IF
+ END IF
+ END IF
+ END IF
+* .. end of the "R**T or R" branch
+ END IF
+*
+* The Q matrix from the first QRF is built into the left singular
+* vectors matrix U.
+*
+ IF ( .NOT. WNTUF )
+ $ CALL SORMQR( 'L', 'N', M, N1, N, A, LDA, WORK, U,
+ $ LDU, WORK(N+1), LWORK-N, IERR )
+ IF ( ROWPRM .AND. .NOT.WNTUF )
+ $ CALL SLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 )
+*
+* ... end of the "full SVD" branch
+ END IF
+*
+* Check whether some singular values are returned as zeros, e.g.
+* due to underflow, and update the numerical rank.
+ p = NR
+ DO 4001 q = p, 1, -1
+ IF ( S(q) .GT. ZERO ) GO TO 4002
+ NR = NR - 1
+ 4001 CONTINUE
+ 4002 CONTINUE
+*
+* .. if numerical rank deficiency is detected, the truncated
+* singular values are set to zero.
+ IF ( NR .LT. N ) CALL SLASET( 'G', N-NR,1, ZERO,ZERO, S(NR+1), N )
+* .. undo scaling; this may cause overflow in the largest singular
+* values.
+ IF ( ASCALED )
+ $ CALL SLASCL( 'G',0,0, ONE,SQRT(REAL(M)), NR,1, S, N, IERR )
+ IF ( CONDA ) RWORK(1) = SCONDA
+ RWORK(2) = p - NR
+* .. p-NR is the number of singular values that are computed as
+* exact zeros in SGESVD() applied to the (possibly truncated)
+* full row rank triangular (trapezoidal) factor of A.
+ NUMRANK = NR
+*
+ RETURN
+*
+* End of SGESVDQ
+*
+ END
diff --git a/lapack-netlib/SRC/sgesvj.f b/lapack-netlib/SRC/sgesvj.f
index 7a7901135..fee5aba4a 100644
--- a/lapack-netlib/SRC/sgesvj.f
+++ b/lapack-netlib/SRC/sgesvj.f
@@ -90,13 +90,13 @@
*> JOBV is CHARACTER*1
*> Specifies whether to compute the right singular vectors, that
*> is, the matrix V:
-*> = 'V' : the matrix V is computed and returned in the array V
-*> = 'A' : the Jacobi rotations are applied to the MV-by-N
+*> = 'V': the matrix V is computed and returned in the array V
+*> = 'A': the Jacobi rotations are applied to the MV-by-N
*> array V. In other words, the right singular vector
*> matrix V is not computed explicitly; instead it is
*> applied to an MV-by-N matrix initially stored in the
*> first MV rows of V.
-*> = 'N' : the matrix V is not computed and the array V is not
+*> = 'N': the matrix V is not computed and the array V is not
*> referenced
*> \endverbatim
*>
@@ -118,8 +118,8 @@
*> A is REAL array, dimension (LDA,N)
*> On entry, the M-by-N matrix A.
*> On exit,
-*> If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C':
-*> If INFO .EQ. 0 :
+*> If JOBU = 'U' .OR. JOBU = 'C':
+*> If INFO = 0:
*> RANKA orthonormal columns of U are returned in the
*> leading RANKA columns of the array A. Here RANKA <= N
*> is the number of computed singular values of A that are
@@ -129,9 +129,9 @@
*> in the array WORK as RANKA=NINT(WORK(2)). Also see the
*> descriptions of SVA and WORK. The computed columns of U
*> are mutually numerically orthogonal up to approximately
-*> TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
+*> TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'),
*> see the description of JOBU.
-*> If INFO .GT. 0,
+*> If INFO > 0,
*> the procedure SGESVJ did not converge in the given number
*> of iterations (sweeps). In that case, the computed
*> columns of U may not be orthogonal up to TOL. The output
@@ -139,8 +139,8 @@
*> values in SVA(1:N)) and V is still a decomposition of the
*> input matrix A in the sense that the residual
*> ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
-*> If JOBU .EQ. 'N':
-*> If INFO .EQ. 0 :
+*> If JOBU = 'N':
+*> If INFO = 0:
*> Note that the left singular vectors are 'for free' in the
*> one-sided Jacobi SVD algorithm. However, if only the
*> singular values are needed, the level of numerical
@@ -149,7 +149,7 @@
*> numerically orthogonal up to approximately M*EPS. Thus,
*> on exit, A contains the columns of U scaled with the
*> corresponding singular values.
-*> If INFO .GT. 0 :
+*> If INFO > 0:
*> the procedure SGESVJ did not converge in the given number
*> of iterations (sweeps).
*> \endverbatim
@@ -164,9 +164,9 @@
*> \verbatim
*> SVA is REAL array, dimension (N)
*> On exit,
-*> If INFO .EQ. 0 :
+*> If INFO = 0 :
*> depending on the value SCALE = WORK(1), we have:
-*> If SCALE .EQ. ONE:
+*> If SCALE = ONE:
*> SVA(1:N) contains the computed singular values of A.
*> During the computation SVA contains the Euclidean column
*> norms of the iterated matrices in the array A.
@@ -175,7 +175,7 @@
*> factored representation is due to the fact that some of the
*> singular values of A might underflow or overflow.
*>
-*> If INFO .GT. 0 :
+*> If INFO > 0 :
*> the procedure SGESVJ did not converge in the given number of
*> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
*> \endverbatim
@@ -183,7 +183,7 @@
*> \param[in] MV
*> \verbatim
*> MV is INTEGER
-*> If JOBV .EQ. 'A', then the product of Jacobi rotations in SGESVJ
+*> If JOBV = 'A', then the product of Jacobi rotations in SGESVJ
*> is applied to the first MV rows of V. See the description of JOBV.
*> \endverbatim
*>
@@ -201,16 +201,16 @@
*> \param[in] LDV
*> \verbatim
*> LDV is INTEGER
-*> The leading dimension of the array V, LDV .GE. 1.
-*> If JOBV .EQ. 'V', then LDV .GE. max(1,N).
-*> If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
+*> The leading dimension of the array V, LDV >= 1.
+*> If JOBV = 'V', then LDV >= max(1,N).
+*> If JOBV = 'A', then LDV >= max(1,MV) .
*> \endverbatim
*>
*> \param[in,out] WORK
*> \verbatim
*> WORK is REAL array, dimension (LWORK)
*> On entry,
-*> If JOBU .EQ. 'C' :
+*> If JOBU = 'C' :
*> WORK(1) = CTOL, where CTOL defines the threshold for convergence.
*> The process stops if all columns of A are mutually
*> orthogonal up to CTOL*EPS, EPS=SLAMCH('E').
@@ -230,7 +230,7 @@
*> WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
*> This is useful information in cases when SGESVJ did
*> not converge, as it can be used to estimate whether
-*> the output is stil useful and for post festum analysis.
+*> the output is still useful and for post festum analysis.
*> WORK(6) = the largest absolute value over all sines of the
*> Jacobi rotation angles in the last sweep. It can be
*> useful for a post festum analysis.
@@ -245,9 +245,9 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0 : successful exit.
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value
-*> > 0 : SGESVJ did not converge in the maximal allowed number (30)
+*> = 0: successful exit.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value
+*> > 0: SGESVJ did not converge in the maximal allowed number (30)
*> of sweeps. The output may still be useful. See the
*> description of WORK.
*> \endverbatim
diff --git a/lapack-netlib/SRC/sgesvxx.f b/lapack-netlib/SRC/sgesvxx.f
index 281f198d5..7cb29d5ab 100644
--- a/lapack-netlib/SRC/sgesvxx.f
+++ b/lapack-netlib/SRC/sgesvxx.f
@@ -411,7 +411,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@@ -447,14 +447,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> NPARAMS is INTEGER
-*> Specifies the number of parameters set in PARAMS. If .LE. 0, the
+*> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> PARAMS is REAL array, dimension NPARAMS
-*> Specifies algorithm parameters. If an entry is .LT. 0.0, then
+*> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.
@@ -462,9 +462,9 @@
*> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
*> refinement or not.
*> Default: 1.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/sgetc2.f b/lapack-netlib/SRC/sgetc2.f
index b0301b953..6bf0a93c6 100644
--- a/lapack-netlib/SRC/sgetc2.f
+++ b/lapack-netlib/SRC/sgetc2.f
@@ -85,7 +85,7 @@
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
-*> > 0: if INFO = k, U(k, k) is likely to produce owerflow if
+*> > 0: if INFO = k, U(k, k) is likely to produce overflow if
*> we try to solve for x in Ax = b. So U is perturbed to
*> avoid the overflow.
*> \endverbatim
diff --git a/lapack-netlib/SRC/sgetsls.f b/lapack-netlib/SRC/sgetsls.f
index 35af66c19..53d2f9431 100644
--- a/lapack-netlib/SRC/sgetsls.f
+++ b/lapack-netlib/SRC/sgetsls.f
@@ -1,3 +1,5 @@
+*> \brief \b SGETSLS
+*
* Definition:
* ===========
*
@@ -154,7 +156,7 @@
*
*> \date June 2017
*
-*> \ingroup doubleGEsolve
+*> \ingroup realGEsolve
*
* =====================================================================
SUBROUTINE SGETSLS( TRANS, M, N, NRHS, A, LDA, B, LDB,
diff --git a/lapack-netlib/SRC/sggesx.f b/lapack-netlib/SRC/sggesx.f
index 3c6273dcf..25691d164 100644
--- a/lapack-netlib/SRC/sggesx.f
+++ b/lapack-netlib/SRC/sggesx.f
@@ -131,10 +131,10 @@
*> \verbatim
*> SENSE is CHARACTER*1
*> Determines which reciprocal condition numbers are computed.
-*> = 'N' : None are computed;
-*> = 'E' : Computed for average of selected eigenvalues only;
-*> = 'V' : Computed for selected deflating subspaces only;
-*> = 'B' : Computed for both.
+*> = 'N': None are computed;
+*> = 'E': Computed for average of selected eigenvalues only;
+*> = 'V': Computed for selected deflating subspaces only;
+*> = 'B': Computed for both.
*> If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
*> \endverbatim
*>
diff --git a/lapack-netlib/SRC/sgsvj0.f b/lapack-netlib/SRC/sgsvj0.f
index e580efc30..d9177d818 100644
--- a/lapack-netlib/SRC/sgsvj0.f
+++ b/lapack-netlib/SRC/sgsvj0.f
@@ -117,7 +117,7 @@
*> \param[in] MV
*> \verbatim
*> MV is INTEGER
-*> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
+*> If JOBV = 'A', then MV rows of V are post-multipled by a
*> sequence of Jacobi rotations.
*> If JOBV = 'N', then MV is not referenced.
*> \endverbatim
@@ -125,9 +125,9 @@
*> \param[in,out] V
*> \verbatim
*> V is REAL array, dimension (LDV,N)
-*> If JOBV .EQ. 'V' then N rows of V are post-multipled by a
+*> If JOBV = 'V' then N rows of V are post-multipled by a
*> sequence of Jacobi rotations.
-*> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
+*> If JOBV = 'A' then MV rows of V are post-multipled by a
*> sequence of Jacobi rotations.
*> If JOBV = 'N', then V is not referenced.
*> \endverbatim
@@ -136,8 +136,8 @@
*> \verbatim
*> LDV is INTEGER
*> The leading dimension of the array V, LDV >= 1.
-*> If JOBV = 'V', LDV .GE. N.
-*> If JOBV = 'A', LDV .GE. MV.
+*> If JOBV = 'V', LDV >= N.
+*> If JOBV = 'A', LDV >= MV.
*> \endverbatim
*>
*> \param[in] EPS
@@ -157,7 +157,7 @@
*> TOL is REAL
*> TOL is the threshold for Jacobi rotations. For a pair
*> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
-*> applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
+*> applied only if ABS(COS(angle(A(:,p),A(:,q)))) > TOL.
*> \endverbatim
*>
*> \param[in] NSWEEP
@@ -175,14 +175,14 @@
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
-*> LWORK is the dimension of WORK. LWORK .GE. M.
+*> LWORK is the dimension of WORK. LWORK >= M.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0 : successful exit.
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value
+*> = 0: successful exit.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value
*> \endverbatim
*
* Authors:
@@ -1045,7 +1045,7 @@
1993 CONTINUE
* end i=1:NSWEEP loop
-* #:) Reaching this point means that the procedure has comleted the given
+* #:) Reaching this point means that the procedure has completed the given
* number of iterations.
INFO = NSWEEP - 1
GO TO 1995
diff --git a/lapack-netlib/SRC/sgsvj1.f b/lapack-netlib/SRC/sgsvj1.f
index 49b81cf4f..ea4ba2e0e 100644
--- a/lapack-netlib/SRC/sgsvj1.f
+++ b/lapack-netlib/SRC/sgsvj1.f
@@ -61,7 +61,7 @@
*> In terms of the columns of A, the first N1 columns are rotated 'against'
*> the remaining N-N1 columns, trying to increase the angle between the
*> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
-*> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
+*> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parameter.
*> The number of sweeps is given in NSWEEP and the orthogonality threshold
*> is given in TOL.
*> \endverbatim
@@ -147,7 +147,7 @@
*> \param[in] MV
*> \verbatim
*> MV is INTEGER
-*> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
+*> If JOBV = 'A', then MV rows of V are post-multipled by a
*> sequence of Jacobi rotations.
*> If JOBV = 'N', then MV is not referenced.
*> \endverbatim
@@ -155,9 +155,9 @@
*> \param[in,out] V
*> \verbatim
*> V is REAL array, dimension (LDV,N)
-*> If JOBV .EQ. 'V' then N rows of V are post-multipled by a
+*> If JOBV = 'V' then N rows of V are post-multipled by a
*> sequence of Jacobi rotations.
-*> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
+*> If JOBV = 'A' then MV rows of V are post-multipled by a
*> sequence of Jacobi rotations.
*> If JOBV = 'N', then V is not referenced.
*> \endverbatim
@@ -166,8 +166,8 @@
*> \verbatim
*> LDV is INTEGER
*> The leading dimension of the array V, LDV >= 1.
-*> If JOBV = 'V', LDV .GE. N.
-*> If JOBV = 'A', LDV .GE. MV.
+*> If JOBV = 'V', LDV >= N.
+*> If JOBV = 'A', LDV >= MV.
*> \endverbatim
*>
*> \param[in] EPS
@@ -187,7 +187,7 @@
*> TOL is REAL
*> TOL is the threshold for Jacobi rotations. For a pair
*> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
-*> applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
+*> applied only if ABS(COS(angle(A(:,p),A(:,q)))) > TOL.
*> \endverbatim
*>
*> \param[in] NSWEEP
@@ -205,14 +205,14 @@
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
-*> LWORK is the dimension of WORK. LWORK .GE. M.
+*> LWORK is the dimension of WORK. LWORK >= M.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0 : successful exit.
-*> < 0 : if INFO = -i, then the i-th argument had an illegal value
+*> = 0: successful exit.
+*> < 0: if INFO = -i, then the i-th argument had an illegal value
*> \endverbatim
*
* Authors:
diff --git a/lapack-netlib/SRC/shseqr.f b/lapack-netlib/SRC/shseqr.f
index 5654a4682..b5707f2c3 100644
--- a/lapack-netlib/SRC/shseqr.f
+++ b/lapack-netlib/SRC/shseqr.f
@@ -70,7 +70,7 @@
*> \param[in] N
*> \verbatim
*> N is INTEGER
-*> The order of the matrix H. N .GE. 0.
+*> The order of the matrix H. N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
@@ -87,7 +87,7 @@
*> set by a previous call to SGEBAL, and then passed to ZGEHRD
*> when the matrix output by SGEBAL is reduced to Hessenberg
*> form. Otherwise ILO and IHI should be set to 1 and N
-*> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
+*> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
*> If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim
*>
@@ -100,20 +100,20 @@
*> (the Schur form); 2-by-2 diagonal blocks (corresponding to
*> complex conjugate pairs of eigenvalues) are returned in
*> standard form, with H(i,i) = H(i+1,i+1) and
-*> H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
+*> H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and JOB = 'E', the
*> contents of H are unspecified on exit. (The output value of
-*> H when INFO.GT.0 is given under the description of INFO
+*> H when INFO > 0 is given under the description of INFO
*> below.)
*>
*> Unlike earlier versions of SHSEQR, this subroutine may
-*> explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
+*> explicitly H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1
*> or j = IHI+1, IHI+2, ... N.
*> \endverbatim
*>
*> \param[in] LDH
*> \verbatim
*> LDH is INTEGER
-*> The leading dimension of the array H. LDH .GE. max(1,N).
+*> The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim
*>
*> \param[out] WR
@@ -128,8 +128,8 @@
*> The real and imaginary parts, respectively, of the computed
*> eigenvalues. If two eigenvalues are computed as a complex
*> conjugate pair, they are stored in consecutive elements of
-*> WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
-*> WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
+*> WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and
+*> WI(i+1) < 0. If JOB = 'S', the eigenvalues are stored in
*> the same order as on the diagonal of the Schur form returned
*> in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
*> diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
@@ -148,7 +148,7 @@
*> if INFO = 0, Z contains Q*Z.
*> Normally Q is the orthogonal matrix generated by SORGHR
*> after the call to SGEHRD which formed the Hessenberg matrix
-*> H. (The output value of Z when INFO.GT.0 is given under
+*> H. (The output value of Z when INFO > 0 is given under
*> the description of INFO below.)
*> \endverbatim
*>
@@ -156,7 +156,7 @@
*> \verbatim
*> LDZ is INTEGER
*> The leading dimension of the array Z. if COMPZ = 'I' or
-*> COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1.
+*> COMPZ = 'V', then LDZ >= MAX(1,N). Otherwise, LDZ >= 1.
*> \endverbatim
*>
*> \param[out] WORK
@@ -169,7 +169,7 @@
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
-*> The dimension of the array WORK. LWORK .GE. max(1,N)
+*> The dimension of the array WORK. LWORK >= max(1,N)
*> is sufficient and delivers very good and sometimes
*> optimal performance. However, LWORK as large as 11*N
*> may be required for optimal performance. A workspace
@@ -187,21 +187,21 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0: successful exit
-*> .LT. 0: if INFO = -i, the i-th argument had an illegal
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal
*> value
-*> .GT. 0: if INFO = i, SHSEQR failed to compute all of
+*> > 0: if INFO = i, SHSEQR failed to compute all of
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
*> and WI contain those eigenvalues which have been
*> successfully computed. (Failures are rare.)
*>
-*> If INFO .GT. 0 and JOB = 'E', then on exit, the
+*> If INFO > 0 and JOB = 'E', then on exit, the
*> remaining unconverged eigenvalues are the eigen-
*> values of the upper Hessenberg matrix rows and
*> columns ILO through INFO of the final, output
*> value of H.
*>
-*> If INFO .GT. 0 and JOB = 'S', then on exit
+*> If INFO > 0 and JOB = 'S', then on exit
*>
*> (*) (initial value of H)*U = U*(final value of H)
*>
@@ -209,19 +209,19 @@
*> value of H is upper Hessenberg and quasi-triangular
*> in rows and columns INFO+1 through IHI.
*>
-*> If INFO .GT. 0 and COMPZ = 'V', then on exit
+*> If INFO > 0 and COMPZ = 'V', then on exit
*>
*> (final value of Z) = (initial value of Z)*U
*>
*> where U is the orthogonal matrix in (*) (regard-
*> less of the value of JOB.)
*>
-*> If INFO .GT. 0 and COMPZ = 'I', then on exit
+*> If INFO > 0 and COMPZ = 'I', then on exit
*> (final value of Z) = U
*> where U is the orthogonal matrix in (*) (regard-
*> less of the value of JOB.)
*>
-*> If INFO .GT. 0 and COMPZ = 'N', then Z is not
+*> If INFO > 0 and COMPZ = 'N', then Z is not
*> accessed.
*> \endverbatim
*
@@ -261,8 +261,8 @@
*> This depends on ILO, IHI and NS. NS is the
*> number of simultaneous shifts returned
*> by ILAENV(ISPEC=15). (See ISPEC=15 below.)
-*> The default for (IHI-ILO+1).LE.500 is NS.
-*> The default for (IHI-ILO+1).GT.500 is 3*NS/2.
+*> The default for (IHI-ILO+1) <= 500 is NS.
+*> The default for (IHI-ILO+1) > 500 is 3*NS/2.
*>
*> ISPEC=14: Nibble crossover point. (See IPARMQ for
*> details.) Default: 14% of deflation window
@@ -341,8 +341,8 @@
PARAMETER ( NTINY = 11 )
*
* ==== NL allocates some local workspace to help small matrices
-* . through a rare SLAHQR failure. NL .GT. NTINY = 11 is
-* . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
+* . through a rare SLAHQR 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/sla_gbrcond.f b/lapack-netlib/SRC/sla_gbrcond.f
index 36aa93dc9..7f2c4062e 100644
--- a/lapack-netlib/SRC/sla_gbrcond.f
+++ b/lapack-netlib/SRC/sla_gbrcond.f
@@ -140,13 +140,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is REAL array, dimension (5*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] IWORK
+*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/sla_gbrfsx_extended.f b/lapack-netlib/SRC/sla_gbrfsx_extended.f
index a81feb45e..0fd1fd350 100644
--- a/lapack-netlib/SRC/sla_gbrfsx_extended.f
+++ b/lapack-netlib/SRC/sla_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/sla_gercond.f b/lapack-netlib/SRC/sla_gercond.f
index 349a1b5be..e54e0d7b4 100644
--- a/lapack-netlib/SRC/sla_gercond.f
+++ b/lapack-netlib/SRC/sla_gercond.f
@@ -122,13 +122,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is REAL array, dimension (3*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] IWORK
+*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N).
*> Workspace.2
diff --git a/lapack-netlib/SRC/sla_gerfsx_extended.f b/lapack-netlib/SRC/sla_gerfsx_extended.f
index 1795ea975..84d1ae31b 100644
--- a/lapack-netlib/SRC/sla_gerfsx_extended.f
+++ b/lapack-netlib/SRC/sla_gerfsx_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
@@ -257,7 +257,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/sla_porcond.f b/lapack-netlib/SRC/sla_porcond.f
index 9dd7c587b..729581f46 100644
--- a/lapack-netlib/SRC/sla_porcond.f
+++ b/lapack-netlib/SRC/sla_porcond.f
@@ -112,13 +112,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is REAL array, dimension (3*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] IWORK
+*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/sla_porfsx_extended.f b/lapack-netlib/SRC/sla_porfsx_extended.f
index 27baa20d1..abbfebb83 100644
--- a/lapack-netlib/SRC/sla_porfsx_extended.f
+++ b/lapack-netlib/SRC/sla_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/sla_syrcond.f b/lapack-netlib/SRC/sla_syrcond.f
index c4b204cc6..0c9e2b361 100644
--- a/lapack-netlib/SRC/sla_syrcond.f
+++ b/lapack-netlib/SRC/sla_syrcond.f
@@ -118,13 +118,13 @@
*> i > 0: The ith argument is invalid.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is REAL array, dimension (3*N).
*> Workspace.
*> \endverbatim
*>
-*> \param[in] IWORK
+*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N).
*> Workspace.
diff --git a/lapack-netlib/SRC/sla_syrfsx_extended.f b/lapack-netlib/SRC/sla_syrfsx_extended.f
index f7b909ac0..a83a9db98 100644
--- a/lapack-netlib/SRC/sla_syrfsx_extended.f
+++ b/lapack-netlib/SRC/sla_syrfsx_extended.f
@@ -67,11 +67,11 @@
*> \verbatim
*> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement.
-*> The value is defined by ILAPREC(P) where P is a CHARACTER and
-*> P = 'S': Single
+*> The value is defined by ILAPREC(P) where P is a CHARACTER and P
+*> = 'S': Single
*> = 'D': Double
*> = 'I': Indigenous
-*> = 'X', 'E': Extra
+*> = 'X' or 'E': Extra
*> \endverbatim
*>
*> \param[in] UPLO
@@ -255,7 +255,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
-*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
+*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
diff --git a/lapack-netlib/SRC/sla_syrpvgrw.f b/lapack-netlib/SRC/sla_syrpvgrw.f
index f5eb81b1f..a0a235ee3 100644
--- a/lapack-netlib/SRC/sla_syrpvgrw.f
+++ b/lapack-netlib/SRC/sla_syrpvgrw.f
@@ -101,7 +101,7 @@
*> as determined by SSYTRF.
*> \endverbatim
*>
-*> \param[in] WORK
+*> \param[out] WORK
*> \verbatim
*> WORK is REAL array, dimension (2*N)
*> \endverbatim
diff --git a/lapack-netlib/SRC/sla_wwaddw.f b/lapack-netlib/SRC/sla_wwaddw.f
index 96a7d3542..e390c9fab 100644
--- a/lapack-netlib/SRC/sla_wwaddw.f
+++ b/lapack-netlib/SRC/sla_wwaddw.f
@@ -36,7 +36,7 @@
*> SLA_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/slaed4.f b/lapack-netlib/SRC/slaed4.f
index c65cba75a..64260843f 100644
--- a/lapack-netlib/SRC/slaed4.f
+++ b/lapack-netlib/SRC/slaed4.f
@@ -82,7 +82,7 @@
*> \param[out] DELTA
*> \verbatim
*> DELTA is REAL array, dimension (N)
-*> If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th
+*> If N > 2, DELTA contains (D(j) - lambda_I) in its j-th
*> component. If N = 1, then DELTA(1) = 1. If N = 2, see SLAED5
*> for detail. The vector DELTA contains the information necessary
*> to construct the eigenvectors by SLAED3 and SLAED9.
diff --git a/lapack-netlib/SRC/slaed8.f b/lapack-netlib/SRC/slaed8.f
index 5ec117cb5..2e3f6f51f 100644
--- a/lapack-netlib/SRC/slaed8.f
+++ b/lapack-netlib/SRC/slaed8.f
@@ -353,7 +353,7 @@
Z( I ) = W( INDX( I ) )
40 CONTINUE
*
-* Calculate the allowable deflation tolerence
+* Calculate the allowable deflation tolerance
*
IMAX = ISAMAX( N, Z, 1 )
JMAX = ISAMAX( N, D, 1 )
diff --git a/lapack-netlib/SRC/slagtf.f b/lapack-netlib/SRC/slagtf.f
index d3f0b6813..59ef097a7 100644
--- a/lapack-netlib/SRC/slagtf.f
+++ b/lapack-netlib/SRC/slagtf.f
@@ -125,7 +125,7 @@
*> then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
*> returns the smallest positive integer j such that
*>
-*> abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
+*> abs( u(j,j) ) <= norm( (T - lambda*I)(j) )*TOL,
*>
*> where norm( A(j) ) denotes the sum of the absolute values of
*> the jth row of the matrix A. If no such j exists then IN(n)
@@ -137,8 +137,8 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0 : successful exit
-*> .lt. 0: if INFO = -k, the kth argument had an illegal value
+*> = 0: successful exit
+*> < 0: if INFO = -k, the kth argument had an illegal value
*> \endverbatim
*
* Authors:
diff --git a/lapack-netlib/SRC/slagts.f b/lapack-netlib/SRC/slagts.f
index 0c3c5239f..e0c8892d7 100644
--- a/lapack-netlib/SRC/slagts.f
+++ b/lapack-netlib/SRC/slagts.f
@@ -122,12 +122,12 @@
*> \param[in,out] TOL
*> \verbatim
*> TOL is REAL
-*> On entry, with JOB .lt. 0, TOL should be the minimum
+*> On entry, with JOB < 0, TOL should be the minimum
*> perturbation to be made to very small diagonal elements of U.
*> TOL should normally be chosen as about eps*norm(U), where eps
*> is the relative machine precision, but if TOL is supplied as
*> non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
-*> If JOB .gt. 0 then TOL is not referenced.
+*> If JOB > 0 then TOL is not referenced.
*>
*> On exit, TOL is changed as described above, only if TOL is
*> non-positive on entry. Otherwise TOL is unchanged.
@@ -136,14 +136,14 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0 : successful exit
-*> .lt. 0: if INFO = -i, the i-th argument had an illegal value
-*> .gt. 0: overflow would occur when computing the INFO(th)
-*> element of the solution vector x. This can only occur
-*> when JOB is supplied as positive and either means
-*> that a diagonal element of U is very small, or that
-*> the elements of the right-hand side vector y are very
-*> large.
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> > 0: overflow would occur when computing the INFO(th)
+*> element of the solution vector x. This can only occur
+*> when JOB is supplied as positive and either means
+*> that a diagonal element of U is very small, or that
+*> the elements of the right-hand side vector y are very
+*> large.
*> \endverbatim
*
* Authors:
diff --git a/lapack-netlib/SRC/slahqr.f b/lapack-netlib/SRC/slahqr.f
index d91826e61..e5642d2bf 100644
--- a/lapack-netlib/SRC/slahqr.f
+++ b/lapack-netlib/SRC/slahqr.f
@@ -150,26 +150,26 @@
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
-*> = 0: successful exit
-*> .GT. 0: If INFO = i, SLAHQR failed to compute all the
+*> = 0: successful exit
+*> > 0: If INFO = i, SLAHQR failed to compute all the
*> eigenvalues ILO to IHI in a total of 30 iterations
*> per eigenvalue; elements i+1:ihi of WR and WI
*> 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, output
+*> 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/slaln2.f b/lapack-netlib/SRC/slaln2.f
index f9ceee7b7..4c6a55ec7 100644
--- a/lapack-netlib/SRC/slaln2.f
+++ b/lapack-netlib/SRC/slaln2.f
@@ -49,7 +49,7 @@
*> the first column of each being the real part and the second
*> being the imaginary part.
*>
-*> "s" is a scaling factor (.LE. 1), computed by SLALN2, which is
+*> "s" is a scaling factor (<= 1), computed by SLALN2, which is
*> so chosen that X can be computed without overflow. X is further
*> scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
*> than overflow.
diff --git a/lapack-netlib/SRC/slamswlq.f b/lapack-netlib/SRC/slamswlq.f
index b13d02b6c..59ab1a6ee 100644
--- a/lapack-netlib/SRC/slamswlq.f
+++ b/lapack-netlib/SRC/slamswlq.f
@@ -1,3 +1,4 @@
+*> \brief \b SLAMSWLQ
*
* Definition:
* ===========
diff --git a/lapack-netlib/SRC/slamtsqr.f b/lapack-netlib/SRC/slamtsqr.f
index 84ac86ee2..58905ab46 100644
--- a/lapack-netlib/SRC/slamtsqr.f
+++ b/lapack-netlib/SRC/slamtsqr.f
@@ -1,3 +1,4 @@
+*> \brief \b SLAMTSQR
*
* Definition:
* ===========
diff --git a/lapack-netlib/SRC/slangb.f b/lapack-netlib/SRC/slangb.f
index fd538b1b7..706e07501 100644
--- a/lapack-netlib/SRC/slangb.f
+++ b/lapack-netlib/SRC/slangb.f
@@ -129,6 +129,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
@@ -139,22 +140,24 @@
*
* =====================================================================
*
-*
* .. Parameters ..
REAL ONE, ZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
* ..
* .. Local Scalars ..
INTEGER I, J, K, L
- REAL SCALE, SUM, VALUE, TEMP
+ REAL SUM, VALUE, TEMP
* ..
-* .. External Subroutines ..
- EXTERNAL SLASSQ
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
EXTERNAL LSAME, SISNAN
* ..
+* .. External Subroutines ..
+ EXTERNAL SLASSQ, SCOMBSSQ
+* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
* ..
@@ -206,15 +209,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 SLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
SLANGB = VALUE
diff --git a/lapack-netlib/SRC/slange.f b/lapack-netlib/SRC/slange.f
index 2eb8d7d14..0c80f1d40 100644
--- a/lapack-netlib/SRC/slange.f
+++ b/lapack-netlib/SRC/slange.f
@@ -119,6 +119,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM
INTEGER LDA, M, N
@@ -135,10 +136,13 @@
* ..
* .. Local Scalars ..
INTEGER I, J
- REAL SCALE, SUM, VALUE, TEMP
+ REAL SUM, VALUE, TEMP
+* ..
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Subroutines ..
- EXTERNAL SLASSQ
+ EXTERNAL SLASSQ, SCOMBSSQ
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
@@ -194,13 +198,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 SLASSQ( M, A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( M, A( 1, J ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
SLANGE = VALUE
diff --git a/lapack-netlib/SRC/slanhs.f b/lapack-netlib/SRC/slanhs.f
index c5a077fbf..8913031a2 100644
--- a/lapack-netlib/SRC/slanhs.f
+++ b/lapack-netlib/SRC/slanhs.f
@@ -113,6 +113,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM
INTEGER LDA, N
@@ -129,15 +130,18 @@
* ..
* .. Local Scalars ..
INTEGER I, J
- REAL SCALE, SUM, VALUE
+ REAL SUM, VALUE
* ..
-* .. External Subroutines ..
- EXTERNAL SLASSQ
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
EXTERNAL LSAME, SISNAN
* ..
+* .. External Subroutines ..
+ EXTERNAL SLASSQ, SCOMBSSQ
+* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT
* ..
@@ -188,13 +192,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 SLASSQ( MIN( N, J+1 ), A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( MIN( N, J+1 ), A( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
SLANHS = VALUE
diff --git a/lapack-netlib/SRC/slansb.f b/lapack-netlib/SRC/slansb.f
index 8f3fe1eb9..23519025d 100644
--- a/lapack-netlib/SRC/slansb.f
+++ b/lapack-netlib/SRC/slansb.f
@@ -134,6 +134,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
@@ -150,15 +151,18 @@
* ..
* .. Local Scalars ..
INTEGER I, J, L
- REAL ABSA, SCALE, SUM, VALUE
+ REAL ABSA, SUM, VALUE
* ..
-* .. External Subroutines ..
- EXTERNAL SLASSQ
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
EXTERNAL LSAME, SISNAN
* ..
+* .. External Subroutines ..
+ EXTERNAL SLASSQ, SCOMBSSQ
+* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
* ..
@@ -225,29 +229,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 SLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ),
- $ 1, SCALE, SUM )
+ $ 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE
L = K + 1
ELSE
DO 120 J = 1, N - 1
- CALL SLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( MIN( N-J, K ), AB( 2, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE
L = 1
END IF
- SUM = 2*SUM
+ SSQ( 2 ) = 2*SSQ( 2 )
ELSE
L = 1
END IF
- CALL SLASSQ( N, AB( L, 1 ), LDAB, SCALE, SUM )
- VALUE = SCALE*SQRT( SUM )
+*
+* Sum diagonal
+*
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( N, AB( L, 1 ), LDAB, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
SLANSB = VALUE
diff --git a/lapack-netlib/SRC/slansp.f b/lapack-netlib/SRC/slansp.f
index 35390cd1c..7e29d778b 100644
--- a/lapack-netlib/SRC/slansp.f
+++ b/lapack-netlib/SRC/slansp.f
@@ -119,6 +119,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM, UPLO
INTEGER N
@@ -135,15 +136,18 @@
* ..
* .. Local Scalars ..
INTEGER I, J, K
- REAL ABSA, SCALE, SUM, VALUE
+ REAL ABSA, SUM, VALUE
* ..
-* .. External Subroutines ..
- EXTERNAL SLASSQ
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
EXTERNAL LSAME, SISNAN
* ..
+* .. External Subroutines ..
+ EXTERNAL SLASSQ, SCOMBSSQ
+* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
* ..
@@ -217,31 +221,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 SLASSQ( J-1, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
K = K + J
110 CONTINUE
ELSE
DO 120 J = 1, N - 1
- CALL SLASSQ( N-J, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( 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( AP( K ).NE.ZERO ) THEN
ABSA = ABS( 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
@@ -250,7 +271,8 @@
K = K + N - I + 1
END IF
130 CONTINUE
- VALUE = SCALE*SQRT( SUM )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
SLANSP = VALUE
diff --git a/lapack-netlib/SRC/slansy.f b/lapack-netlib/SRC/slansy.f
index c8400e530..66ff1c5c7 100644
--- a/lapack-netlib/SRC/slansy.f
+++ b/lapack-netlib/SRC/slansy.f
@@ -127,6 +127,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER NORM, UPLO
INTEGER LDA, N
@@ -143,15 +144,18 @@
* ..
* .. Local Scalars ..
INTEGER I, J
- REAL ABSA, SCALE, SUM, VALUE
+ REAL ABSA, SUM, VALUE
* ..
-* .. External Subroutines ..
- EXTERNAL SLASSQ
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
EXTERNAL LSAME, SISNAN
* ..
+* .. External Subroutines ..
+ EXTERNAL SLASSQ, SCOMBSSQ
+* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
* ..
@@ -216,21 +220,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 SLASSQ( J-1, A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( J-1, A( 1, J ), 1, COLSSQ(1), COLSSQ(2) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE
ELSE
DO 120 J = 1, N - 1
- CALL SLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( N-J, A( J+1, J ), 1, COLSSQ(1), COLSSQ(2) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE
END IF
- SUM = 2*SUM
- CALL SLASSQ( N, A, LDA+1, SCALE, SUM )
- VALUE = SCALE*SQRT( SUM )
+ SSQ( 2 ) = 2*SSQ( 2 )
+*
+* Sum diagonal
+*
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( N, A, LDA+1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
SLANSY = VALUE
diff --git a/lapack-netlib/SRC/slantb.f b/lapack-netlib/SRC/slantb.f
index 3588779cb..5b94618e1 100644
--- a/lapack-netlib/SRC/slantb.f
+++ b/lapack-netlib/SRC/slantb.f
@@ -145,6 +145,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
@@ -162,15 +163,18 @@
* .. Local Scalars ..
LOGICAL UDIAG
INTEGER I, J, L
- REAL SCALE, SUM, VALUE
+ REAL SUM, VALUE
* ..
-* .. External Subroutines ..
- EXTERNAL SLASSQ
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
EXTERNAL LSAME, SISNAN
* ..
+* .. External Subroutines ..
+ EXTERNAL SLASSQ, SCOMBSSQ
+* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
* ..
@@ -311,46 +315,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 SLASSQ( 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 SCOMBSSQ( 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 SLASSQ( MIN( J, K+1 ), AB( MAX( K+2-J, 1 ), J ),
- $ 1, SCALE, SUM )
+ $ 1, COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( 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 SLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( MIN( N-J, K ), AB( 2, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
300 CONTINUE
END IF
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 310 J = 1, N
- CALL SLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
310 CONTINUE
END IF
END IF
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
SLANTB = VALUE
diff --git a/lapack-netlib/SRC/slantp.f b/lapack-netlib/SRC/slantp.f
index 1423f5ca3..ab781deac 100644
--- a/lapack-netlib/SRC/slantp.f
+++ b/lapack-netlib/SRC/slantp.f
@@ -129,6 +129,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
+ IMPLICIT NONE
* .. Scalar Arguments ..
CHARACTER DIAG, NORM, UPLO
INTEGER N
@@ -146,15 +147,18 @@
* .. Local Scalars ..
LOGICAL UDIAG
INTEGER I, J, K
- REAL SCALE, SUM, VALUE
+ REAL SUM, VALUE
* ..
-* .. External Subroutines ..
- EXTERNAL SLASSQ
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
EXTERNAL LSAME, SISNAN
* ..
+* .. External Subroutines ..
+ EXTERNAL SLASSQ, SCOMBSSQ
+* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SQRT
* ..
@@ -306,45 +310,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 SLASSQ( J-1, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( J-1, AP( K ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( 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 SLASSQ( J, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( J, AP( K ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( 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 SLASSQ( N-J, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( N-J, AP( K ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( 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 SLASSQ( N-J+1, AP( K ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( N-J+1, AP( K ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( 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
*
SLANTP = VALUE
diff --git a/lapack-netlib/SRC/slantr.f b/lapack-netlib/SRC/slantr.f
index 63b855892..04d29f537 100644
--- a/lapack-netlib/SRC/slantr.f
+++ b/lapack-netlib/SRC/slantr.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 LDA, M, N
@@ -163,15 +164,18 @@
* .. Local Scalars ..
LOGICAL UDIAG
INTEGER I, J
- REAL SCALE, SUM, VALUE
+ REAL SUM, VALUE
* ..
-* .. External Subroutines ..
- EXTERNAL SLASSQ
+* .. Local Arrays ..
+ REAL SSQ( 2 ), COLSSQ( 2 )
* ..
* .. External Functions ..
LOGICAL LSAME, SISNAN
EXTERNAL LSAME, SISNAN
* ..
+* .. External Subroutines ..
+ EXTERNAL SLASSQ, SCOMBSSQ
+* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT
* ..
@@ -281,7 +285,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
@@ -311,38 +315,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 SLASSQ( MIN( M, J-1 ), A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( MIN( M, J-1 ), A( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
290 CONTINUE
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 300 J = 1, N
- CALL SLASSQ( MIN( M, J ), A( 1, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( MIN( M, J ), A( 1, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( 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 SLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, SCALE,
- $ SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( M-J, A( MIN( M, J+1 ), J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
310 CONTINUE
ELSE
- SCALE = ZERO
- SUM = ONE
+ SSQ( 1 ) = ZERO
+ SSQ( 2 ) = ONE
DO 320 J = 1, N
- CALL SLASSQ( M-J+1, A( J, J ), 1, SCALE, SUM )
+ COLSSQ( 1 ) = ZERO
+ COLSSQ( 2 ) = ONE
+ CALL SLASSQ( M-J+1, A( J, J ), 1,
+ $ COLSSQ( 1 ), COLSSQ( 2 ) )
+ CALL SCOMBSSQ( SSQ, COLSSQ )
320 CONTINUE
END IF
END IF
- VALUE = SCALE*SQRT( SUM )
+ VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF
*
SLANTR = VALUE
diff --git a/lapack-netlib/SRC/slanv2.f b/lapack-netlib/SRC/slanv2.f
index e73e5455c..1163446fa 100644
--- a/lapack-netlib/SRC/slanv2.f
+++ b/lapack-netlib/SRC/slanv2.f
@@ -161,7 +161,6 @@
IF( C.EQ.ZERO ) THEN
CS = ONE
SN = ZERO
- GO TO 10
*
ELSE IF( B.EQ.ZERO ) THEN
*
@@ -174,12 +173,12 @@
A = TEMP
B = -C
C = ZERO
- GO TO 10
+*
ELSE IF( (A-D).EQ.ZERO .AND. SIGN( ONE, B ).NE.
$ SIGN( ONE, C ) ) THEN
CS = ONE
SN = ZERO
- GO TO 10
+*
ELSE
*
TEMP = A - D
@@ -207,6 +206,7 @@
SN = C / TAU
B = B - C
C = ZERO
+*
ELSE
*
* Complex eigenvalues, or real (almost) equal eigenvalues.
@@ -268,8 +268,6 @@
END IF
*
END IF
-*
- 10 CONTINUE
*
* Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
*
diff --git a/lapack-netlib/SRC/slaorhr_col_getrfnp.f b/lapack-netlib/SRC/slaorhr_col_getrfnp.f
new file mode 100644
index 000000000..6cc59e538
--- /dev/null
+++ b/lapack-netlib/SRC/slaorhr_col_getrfnp.f
@@ -0,0 +1,248 @@
+*> \brief \b SLAORHR_COL_GETRFNP
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download SLAORHR_COL_GETRFNP + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE SLAORHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
+*
+* .. Scalar Arguments ..
+* INTEGER INFO, LDA, M, N
+* ..
+* .. Array Arguments ..
+* REAL A( LDA, * ), D( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> SLAORHR_COL_GETRFNP computes the modified LU factorization without
+*> pivoting of a real 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 SORHR_COL. In SORHR_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 SLAORHR_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 REAL 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 REAL 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 plus or minus one.
+*> \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 realGEcomputational
+*
+*> \par Contributors:
+* ==================
+*>
+*> \verbatim
+*>
+*> November 2019, Igor Kozachenko,
+*> Computer Science Division,
+*> University of California, Berkeley
+*>
+*> \endverbatim
+*
+* =====================================================================
+ SUBROUTINE SLAORHR_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 ..
+ REAL A( LDA, * ), D( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ONE
+ PARAMETER ( ONE = 1.0E+0 )
+* ..
+* .. Local Scalars ..
+ INTEGER IINFO, J, JB, NB
+* ..
+* .. External Subroutines ..
+ EXTERNAL SGEMM, SLAORHR_COL_GETRFNP2, STRSM, 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( 'SLAORHR_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, 'SLAORHR_COL_GETRFNP', ' ', M, N, -1, -1 )
+
+ IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
+*
+* Use unblocked code.
+*
+ CALL SLAORHR_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 SLAORHR_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 STRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
+ $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
+ $ LDA )
+ IF( J+JB.LE.M ) THEN
+*
+* Update trailing submatrix.
+*
+ CALL SGEMM( 'No transpose', 'No transpose', M-J-JB+1,
+ $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
+ $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
+ $ LDA )
+ END IF
+ END IF
+ END DO
+ END IF
+ RETURN
+*
+* End of SLAORHR_COL_GETRFNP
+*
+ END
\ No newline at end of file
diff --git a/lapack-netlib/SRC/slaorhr_col_getrfnp2.f b/lapack-netlib/SRC/slaorhr_col_getrfnp2.f
new file mode 100644
index 000000000..de604602f
--- /dev/null
+++ b/lapack-netlib/SRC/slaorhr_col_getrfnp2.f
@@ -0,0 +1,305 @@
+*> \brief \b SLAORHR_COL_GETRFNP2
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DLAORHR_GETRF2NP + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* RECURSIVE SUBROUTINE SLAORHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
+*
+* .. Scalar Arguments ..
+* INTEGER INFO, LDA, M, N
+* ..
+* .. Array Arguments ..
+* REAL A( LDA, * ), D( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> SLAORHR_COL_GETRFNP2 computes the modified LU factorization without
+*> pivoting of a real 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 SORHR_COL. In SORHR_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].
+*>
+*> SLAORHR_COL_GETRFNP2 is called to factorize a block by the blocked
+*> routine SLAORHR_COL_GETRFNP, which uses blocked code calling
+*. Level 3 BLAS to update the submatrix. However, SLAORHR_COL_GETRFNP2
+*> is self-sufficient and can be used without SLAORHR_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 REAL 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 REAL 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 plus or minus one.
+*> \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 realGEcomputational
+*
+*> \par Contributors:
+* ==================
+*>
+*> \verbatim
+*>
+*> November 2019, Igor Kozachenko,
+*> Computer Science Division,
+*> University of California, Berkeley
+*>
+*> \endverbatim
+*
+* =====================================================================
+ RECURSIVE SUBROUTINE SLAORHR_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 ..
+ REAL A( LDA, * ), D( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ONE
+ PARAMETER ( ONE = 1.0E+0 )
+* ..
+* .. Local Scalars ..
+ REAL SFMIN
+ INTEGER I, IINFO, N1, N2
+* ..
+* .. External Functions ..
+ REAL SLAMCH
+ EXTERNAL SLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL SGEMM, SSCAL, STRSM, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, SIGN, 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( 'SLAORHR_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 ) = -SIGN( ONE, 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 ) = -SIGN( ONE, 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 = SLAMCH('S')
+*
+* Construct the subdiagonal elements of L
+*
+ IF( ABS( A( 1, 1 ) ) .GE. SFMIN ) THEN
+ CALL SSCAL( M-1, ONE / 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 SLAORHR_COL_GETRFNP2( N1, N1, A, LDA, D, IINFO )
+*
+* Solve for B21
+*
+ CALL STRSM( 'R', 'U', 'N', 'N', M-N1, N1, ONE, A, LDA,
+ $ A( N1+1, 1 ), LDA )
+*
+* Solve for B12
+*
+ CALL STRSM( 'L', 'L', 'N', 'U', N1, N2, ONE, A, LDA,
+ $ A( 1, N1+1 ), LDA )
+*
+* Update B22, i.e. compute the Schur complement
+* B22 := B22 - B21*B12
+*
+ CALL SGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( N1+1, 1 ), LDA,
+ $ A( 1, N1+1 ), LDA, ONE, A( N1+1, N1+1 ), LDA )
+*
+* Factor B22, recursive call
+*
+ CALL SLAORHR_COL_GETRFNP2( M-N1, N2, A( N1+1, N1+1 ), LDA,
+ $ D( N1+1 ), IINFO )
+*
+ END IF
+ RETURN
+*
+* End of SLAORHR_COL_GETRFNP2
+*
+ END