LAPACK: SRC/sgejsv.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
00002      &                   M, N, A, LDA, SVA, U, LDU, V, LDV,
00003      &                   WORK, LWORK, IWORK, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.3.0)                                    --
00006 *
00007 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
00008 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
00009 *     November 2010
00010 *
00011 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00012 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00013 *
00014 * This routine is also part of SIGMA (version 1.23, October 23. 2008.)
00015 * SIGMA is a library of algorithms for highly accurate algorithms for
00016 * computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the
00017 * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0.
00018 *
00019 *     .. Scalar Arguments ..
00020 *
00021       IMPLICIT    NONE
00022       INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
00023 *
00024 *     .. Array Arguments ..
00025 *
00026       REAL        A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
00027      &            WORK( LWORK )
00028       INTEGER     IWORK( * )
00029       CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
00030 *     ..
00031 *
00032 *  Purpose
00033 *  =======
00034 *  SGEJSV computes the singular value decomposition (SVD) of a real M-by-N
00035 *  matrix [A], where M >= N. The SVD of [A] is written as
00036 *
00037 *               [A] = [U] * [SIGMA] * [V]^t,
00038 *
00039 *  where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
00040 *  diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and
00041 *  [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are
00042 *  the singular values of [A]. The columns of [U] and [V] are the left and
00043 *  the right singular vectors of [A], respectively. The matrices [U] and [V]
00044 *  are computed and stored in the arrays U and V, respectively. The diagonal
00045 *  of [SIGMA] is computed and stored in the array SVA.
00046 *
00047 *  Arguments
00048 *  =========
00049 *
00050 *  JOBA   (input) CHARACTER*1
00051 *         Specifies the level of accuracy:
00052 *       = 'C': This option works well (high relative accuracy) if A = B * D,
00053 *              with well-conditioned B and arbitrary diagonal matrix D.
00054 *              The accuracy cannot be spoiled by COLUMN scaling. The
00055 *              accuracy of the computed output depends on the condition of
00056 *              B, and the procedure aims at the best theoretical accuracy.
00057 *              The relative error max_{i=1:N}|d sigma_i| / sigma_i is
00058 *              bounded by f(M,N)*epsilon* cond(B), independent of D.
00059 *              The input matrix is preprocessed with the QRF with column
00060 *              pivoting. This initial preprocessing and preconditioning by
00061 *              a rank revealing QR factorization is common for all values of
00062 *              JOBA. Additional actions are specified as follows:
00063 *       = 'E': Computation as with 'C' with an additional estimate of the
00064 *              condition number of B. It provides a realistic error bound.
00065 *       = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
00066 *              D1, D2, and well-conditioned matrix C, this option gives
00067 *              higher accuracy than the 'C' option. If the structure of the
00068 *              input matrix is not known, and relative accuracy is
00069 *              desirable, then this option is advisable. The input matrix A
00070 *              is preprocessed with QR factorization with FULL (row and
00071 *              column) pivoting.
00072 *       = 'G'  Computation as with 'F' with an additional estimate of the
00073 *              condition number of B, where A=D*B. If A has heavily weighted
00074 *              rows, then using this condition number gives too pessimistic
00075 *              error bound.
00076 *       = 'A': Small singular values are the noise and the matrix is treated
00077 *              as numerically rank defficient. The error in the computed
00078 *              singular values is bounded by f(m,n)*epsilon*||A||.
00079 *              The computed SVD A = U * S * V^t restores A up to
00080 *              f(m,n)*epsilon*||A||.
00081 *              This gives the procedure the licence to discard (set to zero)
00082 *              all singular values below N*epsilon*||A||.
00083 *       = 'R': Similar as in 'A'. Rank revealing property of the initial
00084 *              QR factorization is used do reveal (using triangular factor)
00085 *              a gap sigma_{r+1} < epsilon * sigma_r in which case the
00086 *              numerical RANK is declared to be r. The SVD is computed with
00087 *              absolute error bounds, but more accurately than with 'A'.
00088 * 
00089 *  JOBU   (input) CHARACTER*1
00090 *         Specifies whether to compute the columns of U:
00091 *       = 'U': N columns of U are returned in the array U.
00092 *       = 'F': full set of M left sing. vectors is returned in the array U.
00093 *       = 'W': U may be used as workspace of length M*N. See the description
00094 *              of U.
00095 *       = 'N': U is not computed.
00096 * 
00097 *  JOBV   (input) CHARACTER*1
00098 *         Specifies whether to compute the matrix V:
00099 *       = 'V': N columns of V are returned in the array V; Jacobi rotations
00100 *              are not explicitly accumulated.
00101 *       = 'J': N columns of V are returned in the array V, but they are
00102 *              computed as the product of Jacobi rotations. This option is
00103 *              allowed only if JOBU .NE. 'N', i.e. in computing the full SVD.
00104 *       = 'W': V may be used as workspace of length N*N. See the description
00105 *              of V.
00106 *       = 'N': V is not computed.
00107 * 
00108 *  JOBR   (input) CHARACTER*1
00109 *         Specifies the RANGE for the singular values. Issues the licence to
00110 *         set to zero small positive singular values if they are outside
00111 *         specified range. If A .NE. 0 is scaled so that the largest singular
00112 *         value of c*A is around SQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
00113 *         the licence to kill columns of A whose norm in c*A is less than
00114 *         SQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
00115 *         where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
00116 *       = 'N': Do not kill small columns of c*A. This option assumes that
00117 *              BLAS and QR factorizations and triangular solvers are
00118 *              implemented to work in that range. If the condition of A
00119 *              is greater than BIG, use SGESVJ.
00120 *       = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)]
00121 *              (roughly, as described above). This option is recommended.
00122 *                                             ===========================
00123 *         For computing the singular values in the FULL range [SFMIN,BIG]
00124 *         use SGESVJ.
00125 * 
00126 *  JOBT   (input) CHARACTER*1
00127 *         If the matrix is square then the procedure may determine to use
00128 *         transposed A if A^t seems to be better with respect to convergence.
00129 *         If the matrix is not square, JOBT is ignored. This is subject to
00130 *         changes in the future.
00131 *         The decision is based on two values of entropy over the adjoint
00132 *         orbit of A^t * A. See the descriptions of WORK(6) and WORK(7).
00133 *       = 'T': transpose if entropy test indicates possibly faster
00134 *         convergence of Jacobi process if A^t is taken as input. If A is
00135 *         replaced with A^t, then the row pivoting is included automatically.
00136 *       = 'N': do not speculate.
00137 *         This option can be used to compute only the singular values, or the
00138 *         full SVD (U, SIGMA and V). For only one set of singular vectors
00139 *         (U or V), the caller should provide both U and V, as one of the
00140 *         matrices is used as workspace if the matrix A is transposed.
00141 *         The implementer can easily remove this constraint and make the
00142 *         code more complicated. See the descriptions of U and V.
00143 * 
00144 *  JOBP   (input) CHARACTER*1
00145 *         Issues the licence to introduce structured perturbations to drown
00146 *         denormalized numbers. This licence should be active if the
00147 *         denormals are poorly implemented, causing slow computation,
00148 *         especially in cases of fast convergence (!). For details see [1,2].
00149 *         For the sake of simplicity, this perturbations are included only
00150 *         when the full SVD or only the singular values are requested. The
00151 *         implementer/user can easily add the perturbation for the cases of
00152 *         computing one set of singular vectors.
00153 *       = 'P': introduce perturbation
00154 *       = 'N': do not perturb
00155 *
00156 *  M      (input) INTEGER
00157 *         The number of rows of the input matrix A.  M >= 0.
00158 *
00159 *  N      (input) INTEGER
00160 *         The number of columns of the input matrix A. M >= N >= 0.
00161 *
00162 *  A       (input/workspace) REAL array, dimension (LDA,N)
00163 *          On entry, the M-by-N matrix A.
00164 *
00165 *  LDA     (input) INTEGER
00166 *          The leading dimension of the array A.  LDA >= max(1,M).
00167 *
00168 *  SVA     (workspace/output) REAL array, dimension (N)
00169 *          On exit,
00170 *          - For WORK(1)/WORK(2) = ONE: The singular values of A. During the
00171 *            computation SVA contains Euclidean column norms of the
00172 *            iterated matrices in the array A.
00173 *          - For WORK(1) .NE. WORK(2): The singular values of A are
00174 *            (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if
00175 *            sigma_max(A) overflows or if small singular values have been
00176 *            saved from underflow by scaling the input matrix A.
00177 *          - If JOBR='R' then some of the singular values may be returned
00178 *            as exact zeros obtained by "set to zero" because they are
00179 *            below the numerical rank threshold or are denormalized numbers.
00180 *
00181 *  U       (workspace/output) REAL array, dimension ( LDU, N )
00182 *          If JOBU = 'U', then U contains on exit the M-by-N matrix of
00183 *                         the left singular vectors.
00184 *          If JOBU = 'F', then U contains on exit the M-by-M matrix of
00185 *                         the left singular vectors, including an ONB
00186 *                         of the orthogonal complement of the Range(A).
00187 *          If JOBU = 'W'  .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
00188 *                         then U is used as workspace if the procedure
00189 *                         replaces A with A^t. In that case, [V] is computed
00190 *                         in U as left singular vectors of A^t and then
00191 *                         copied back to the V array. This 'W' option is just
00192 *                         a reminder to the caller that in this case U is
00193 *                         reserved as workspace of length N*N.
00194 *          If JOBU = 'N'  U is not referenced.
00195 *
00196 * LDU      (input) INTEGER
00197 *          The leading dimension of the array U,  LDU >= 1.
00198 *          IF  JOBU = 'U' or 'F' or 'W',  then LDU >= M.
00199 *
00200 *  V       (workspace/output) REAL array, dimension ( LDV, N )
00201 *          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
00202 *                         the right singular vectors;
00203 *          If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
00204 *                         then V is used as workspace if the pprocedure
00205 *                         replaces A with A^t. In that case, [U] is computed
00206 *                         in V as right singular vectors of A^t and then
00207 *                         copied back to the U array. This 'W' option is just
00208 *                         a reminder to the caller that in this case V is
00209 *                         reserved as workspace of length N*N.
00210 *          If JOBV = 'N'  V is not referenced.
00211 *
00212 *  LDV     (input) INTEGER
00213 *          The leading dimension of the array V,  LDV >= 1.
00214 *          If JOBV = 'V' or 'J' or 'W', then LDV >= N.
00215 *
00216 *  WORK    (workspace/output) REAL array, dimension at least LWORK.
00217 *          On exit,
00218 *          WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
00219 *                    that SCALE*SVA(1:N) are the computed singular values
00220 *                    of A. (See the description of SVA().)
00221 *          WORK(2) = See the description of WORK(1).
00222 *          WORK(3) = SCONDA is an estimate for the condition number of
00223 *                    column equilibrated A. (If JOBA .EQ. 'E' or 'G')
00224 *                    SCONDA is an estimate of SQRT(||(R^t * R)^(-1)||_1).
00225 *                    It is computed using SPOCON. It holds
00226 *                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
00227 *                    where R is the triangular factor from the QRF of A.
00228 *                    However, if R is truncated and the numerical rank is
00229 *                    determined to be strictly smaller than N, SCONDA is
00230 *                    returned as -1, thus indicating that the smallest
00231 *                    singular values might be lost.
00232 *
00233 *          If full SVD is needed, the following two condition numbers are
00234 *          useful for the analysis of the algorithm. They are provied for
00235 *          a developer/implementer who is familiar with the details of
00236 *          the method.
00237 *
00238 *          WORK(4) = an estimate of the scaled condition number of the
00239 *                    triangular factor in the first QR factorization.
00240 *          WORK(5) = an estimate of the scaled condition number of the
00241 *                    triangular factor in the second QR factorization.
00242 *          The following two parameters are computed if JOBT .EQ. 'T'.
00243 *          They are provided for a developer/implementer who is familiar
00244 *          with the details of the method.
00245 *
00246 *          WORK(6) = the entropy of A^t*A :: this is the Shannon entropy
00247 *                    of diag(A^t*A) / Trace(A^t*A) taken as point in the
00248 *                    probability simplex.
00249 *          WORK(7) = the entropy of A*A^t.
00250 *
00251 *  LWORK   (input) INTEGER
00252 *          Length of WORK to confirm proper allocation of work space.
00253 *          LWORK depends on the job:
00254 *
00255 *          If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
00256 *            -> .. no scaled condition estimate required ( JOBE.EQ.'N'):
00257 *               LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
00258 *               For optimal performance (blocked code) the optimal value
00259 *               is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
00260 *               block size for xGEQP3/xGEQRF.
00261 *            -> .. an estimate of the scaled condition number of A is
00262 *               required (JOBA='E', 'G'). In this case, LWORK is the maximum
00263 *               of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4N,7).
00264 *
00265 *          If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
00266 *            -> the minimal requirement is LWORK >= max(2*N+M,7).
00267 *            -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7),
00268 *               where NB is the optimal block size.
00269 *
00270 *          If SIGMA and the left singular vectors are needed
00271 *            -> the minimal requirement is LWORK >= max(2*N+M,7).
00272 *            -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7),
00273 *               where NB is the optimal block size.
00274 *
00275 *          If full SVD is needed ( JOBU.EQ.'U' or 'F', JOBV.EQ.'V' ) and
00276 *            -> .. the singular vectors are computed without explicit
00277 *               accumulation of the Jacobi rotations, LWORK >= 6*N+2*N*N
00278 *            -> .. in the iterative part, the Jacobi rotations are
00279 *               explicitly accumulated (option, see the description of JOBV),
00280 *               then the minimal requirement is LWORK >= max(M+3*N+N*N,7).
00281 *               For better performance, if NB is the optimal block size,
00282 *               LWORK >= max(3*N+N*N+M,3*N+N*N+N*NB,7).
00283 *
00284 *  IWORK   (workspace/output) INTEGER array, dimension M+3*N.
00285 *          On exit,
00286 *          IWORK(1) = the numerical rank determined after the initial
00287 *                     QR factorization with pivoting. See the descriptions
00288 *                     of JOBA and JOBR.
00289 *          IWORK(2) = the number of the computed nonzero singular values
00290 *          IWORK(3) = if nonzero, a warning message:
00291 *                     If IWORK(3).EQ.1 then some of the column norms of A
00292 *                     were denormalized floats. The requested high accuracy
00293 *                     is not warranted by the data.
00294 *
00295 *  INFO    (output) INTEGER
00296 *           < 0  : if INFO = -i, then the i-th argument had an illegal value.
00297 *           = 0 :  successfull exit;
00298 *           > 0 :  SGEJSV  did not converge in the maximal allowed number
00299 *                  of sweeps. The computed values may be inaccurate.
00300 *
00301 *  Further Details
00302 *  ===============
00303 *
00304 *  SGEJSV implements a preconditioned Jacobi SVD algorithm. It uses SGEQP3,
00305 *  SGEQRF, and SGELQF as preprocessors and preconditioners. Optionally, an
00306 *  additional row pivoting can be used as a preprocessor, which in some
00307 *  cases results in much higher accuracy. An example is matrix A with the
00308 *  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
00309 *  diagonal matrices and C is well-conditioned matrix. In that case, complete
00310 *  pivoting in the first QR factorizations provides accuracy dependent on the
00311 *  condition number of C, and independent of D1, D2. Such higher accuracy is
00312 *  not completely understood theoretically, but it works well in practice.
00313 *  Further, if A can be written as A = B*D, with well-conditioned B and some
00314 *  diagonal D, then the high accuracy is guaranteed, both theoretically and
00315 *  in software, independent of D. For more details see [1], [2].
00316 *     The computational range for the singular values can be the full range
00317 *  ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
00318 *  & LAPACK routines called by SGEJSV are implemented to work in that range.
00319 *  If that is not the case, then the restriction for safe computation with
00320 *  the singular values in the range of normalized IEEE numbers is that the
00321 *  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
00322 *  overflow. This code (SGEJSV) is best used in this restricted range,
00323 *  meaning that singular values of magnitude below ||A||_2 / SLAMCH('O') are
00324 *  returned as zeros. See JOBR for details on this.
00325 *     Further, this implementation is somewhat slower than the one described
00326 *  in [1,2] due to replacement of some non-LAPACK components, and because
00327 *  the choice of some tuning parameters in the iterative part (SGESVJ) is
00328 *  left to the implementer on a particular machine.
00329 *     The rank revealing QR factorization (in this code: SGEQP3) should be
00330 *  implemented as in [3]. We have a new version of SGEQP3 under development
00331 *  that is more robust than the current one in LAPACK, with a cleaner cut in
00332 *  rank defficient cases. It will be available in the SIGMA library [4].
00333 *  If M is much larger than N, it is obvious that the inital QRF with
00334 *  column pivoting can be preprocessed by the QRF without pivoting. That
00335 *  well known trick is not used in SGEJSV because in some cases heavy row
00336 *  weighting can be treated with complete pivoting. The overhead in cases
00337 *  M much larger than N is then only due to pivoting, but the benefits in
00338 *  terms of accuracy have prevailed. The implementer/user can incorporate
00339 *  this extra QRF step easily. The implementer can also improve data movement
00340 *  (matrix transpose, matrix copy, matrix transposed copy) - this
00341 *  implementation of SGEJSV uses only the simplest, naive data movement.
00342 *
00343 *  Contributors
00344 *
00345 *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
00346 *
00347 *  References
00348 *
00349 * [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
00350 *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
00351 *     LAPACK Working note 169.
00352 * [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
00353 *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
00354 *     LAPACK Working note 170.
00355 * [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
00356 *     factorization software - a case study.
00357 *     ACM Trans. math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
00358 *     LAPACK Working note 176.
00359 * [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
00360 *     QSVD, (H,K)-SVD computations.
00361 *     Department of Mathematics, University of Zagreb, 2008.
00362 *
00363 *  Bugs, examples and comments
00364 *
00365 *  Please report all bugs and send interesting examples and/or comments to
00366 *  drmac@math.hr. Thank you.
00367 *
00368 *  ===========================================================================
00369 *
00370 *     .. Local Parameters ..
00371       REAL        ZERO,         ONE
00372       PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
00373 *     ..
00374 *     .. Local Scalars ..
00375       REAL    AAPP,   AAQQ,   AATMAX, AATMIN, BIG,    BIG1,   COND_OK,
00376      &        CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,  MAXPRJ, SCALEM,
00377      &        SCONDA, SFMIN,  SMALL,  TEMP1,  USCAL1, USCAL2, XSC
00378       INTEGER IERR,   N1,     NR,     NUMRANK,        p, q,   WARNING
00379       LOGICAL ALMORT, DEFR,   ERREST, GOSCAL, JRACC,  KILL,   LSVEC,
00380      &        L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN,
00381      &        NOSCAL, ROWPIV, RSVEC,  TRANSP
00382 *     ..
00383 *     .. Intrinsic Functions ..
00384       INTRINSIC ABS,  ALOG, AMAX1, AMIN1, FLOAT,
00385      &          MAX0, MIN0, NINT,  SIGN,  SQRT
00386 *     ..
00387 *     .. External Functions ..
00388       REAL      SLAMCH, SNRM2
00389       INTEGER   ISAMAX
00390       LOGICAL   LSAME
00391       EXTERNAL  ISAMAX, LSAME, SLAMCH, SNRM2
00392 *     ..
00393 *     .. External Subroutines ..
00394       EXTERNAL  SCOPY,  SGELQF, SGEQP3, SGEQRF, SLACPY, SLASCL,
00395      &          SLASET, SLASSQ, SLASWP, SORGQR, SORMLQ,
00396      &          SORMQR, SPOCON, SSCAL,  SSWAP,  STRSM,  XERBLA
00397 *
00398       EXTERNAL  SGESVJ
00399 *     ..
00400 *
00401 *     Test the input arguments
00402 *
00403       LSVEC  = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
00404       JRACC  = LSAME( JOBV, 'J' )
00405       RSVEC  = LSAME( JOBV, 'V' ) .OR. JRACC
00406       ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
00407       L2RANK = LSAME( JOBA, 'R' )
00408       L2ABER = LSAME( JOBA, 'A' )
00409       ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
00410       L2TRAN = LSAME( JOBT, 'T' )
00411       L2KILL = LSAME( JOBR, 'R' )
00412       DEFR   = LSAME( JOBR, 'N' )
00413       L2PERT = LSAME( JOBP, 'P' )
00414 *
00415       IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
00416      &     ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
00417          INFO = - 1
00418       ELSE IF ( .NOT.( LSVEC  .OR. LSAME( JOBU, 'N' ) .OR.
00419      &                             LSAME( JOBU, 'W' )) ) THEN
00420          INFO = - 2
00421       ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
00422      &   LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN
00423          INFO = - 3
00424       ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) )    THEN
00425          INFO = - 4
00426       ELSE IF ( .NOT. ( L2TRAN .OR. LSAME( JOBT, 'N' ) ) ) THEN
00427          INFO = - 5
00428       ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
00429          INFO = - 6
00430       ELSE IF ( M .LT. 0 ) THEN
00431          INFO = - 7
00432       ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
00433          INFO = - 8
00434       ELSE IF ( LDA .LT. M ) THEN
00435          INFO = - 10
00436       ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
00437          INFO = - 13
00438       ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
00439          INFO = - 14
00440       ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND.
00441      &                           (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR.
00442      & (.NOT.(LSVEC .OR. LSVEC) .AND. ERREST .AND.
00443      &                         (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR.
00444      & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR.
00445      & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR.
00446      & (LSVEC .AND. RSVEC .AND. .NOT.JRACC .AND. (LWORK.LT.6*N+2*N*N))
00447      & .OR. (LSVEC.AND.RSVEC.AND.JRACC.AND.LWORK.LT.MAX0(7,M+3*N+N*N)))
00448      &   THEN
00449          INFO = - 17
00450       ELSE
00451 *        #:)
00452          INFO = 0
00453       END IF
00454 *
00455       IF ( INFO .NE. 0 ) THEN
00456 *       #:(
00457          CALL XERBLA( 'SGEJSV', - INFO )
00458       END IF
00459 *
00460 *     Quick return for void matrix (Y3K safe)
00461 * #:)
00462       IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
00463 *
00464 *     Determine whether the matrix U should be M x N or M x M
00465 *
00466       IF ( LSVEC ) THEN
00467          N1 = N
00468          IF ( LSAME( JOBU, 'F' ) ) N1 = M
00469       END IF
00470 *
00471 *     Set numerical parameters
00472 *
00473 *!    NOTE: Make sure SLAMCH() does not fail on the target architecture.
00474 *
00475       EPSLN = SLAMCH('Epsilon')
00476       SFMIN = SLAMCH('SafeMinimum')
00477       SMALL = SFMIN / EPSLN
00478       BIG   = SLAMCH('O')
00479 *
00480 *     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
00481 *
00482 *(!)  If necessary, scale SVA() to protect the largest norm from
00483 *     overflow. It is possible that this scaling pushes the smallest
00484 *     column norm left from the underflow threshold (extreme case).
00485 *
00486       SCALEM  = ONE / SQRT(FLOAT(M)*FLOAT(N))
00487       NOSCAL  = .TRUE.
00488       GOSCAL  = .TRUE.
00489       DO 1874 p = 1, N
00490          AAPP = ZERO
00491          AAQQ = ONE
00492          CALL SLASSQ( M, A(1,p), 1, AAPP, AAQQ )
00493          IF ( AAPP .GT. BIG ) THEN
00494             INFO = - 9
00495             CALL XERBLA( 'SGEJSV', -INFO )
00496             RETURN
00497          END IF
00498          AAQQ = SQRT(AAQQ)
00499          IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL  ) THEN
00500             SVA(p)  = AAPP * AAQQ
00501          ELSE
00502             NOSCAL  = .FALSE.
00503             SVA(p)  = AAPP * ( AAQQ * SCALEM )
00504             IF ( GOSCAL ) THEN
00505                GOSCAL = .FALSE.
00506                CALL SSCAL( p-1, SCALEM, SVA, 1 )
00507             END IF
00508          END IF
00509  1874 CONTINUE
00510 *
00511       IF ( NOSCAL ) SCALEM = ONE
00512 *
00513       AAPP = ZERO
00514       AAQQ = BIG
00515       DO 4781 p = 1, N
00516          AAPP = AMAX1( AAPP, SVA(p) )
00517          IF ( SVA(p) .NE. ZERO ) AAQQ = AMIN1( AAQQ, SVA(p) )
00518  4781 CONTINUE
00519 *
00520 *     Quick return for zero M x N matrix
00521 * #:)
00522       IF ( AAPP .EQ. ZERO ) THEN
00523          IF ( LSVEC ) CALL SLASET( 'G', M, N1, ZERO, ONE, U, LDU )
00524          IF ( RSVEC ) CALL SLASET( 'G', N, N,  ZERO, ONE, V, LDV )
00525          WORK(1) = ONE
00526          WORK(2) = ONE
00527          IF ( ERREST ) WORK(3) = ONE
00528          IF ( LSVEC .AND. RSVEC ) THEN
00529             WORK(4) = ONE
00530             WORK(5) = ONE
00531          END IF
00532          IF ( L2TRAN ) THEN
00533             WORK(6) = ZERO
00534             WORK(7) = ZERO
00535          END IF
00536          IWORK(1) = 0
00537          IWORK(2) = 0
00538          RETURN
00539       END IF
00540 *
00541 *     Issue warning if denormalized column norms detected. Override the
00542 *     high relative accuracy request. Issue licence to kill columns
00543 *     (set them to zero) whose norm is less than sigma_max / BIG (roughly).
00544 * #:(
00545       WARNING = 0
00546       IF ( AAQQ .LE. SFMIN ) THEN
00547          L2RANK = .TRUE.
00548          L2KILL = .TRUE.
00549          WARNING = 1
00550       END IF
00551 *
00552 *     Quick return for one-column matrix
00553 * #:)
00554       IF ( N .EQ. 1 ) THEN
00555 *
00556          IF ( LSVEC ) THEN
00557             CALL SLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
00558             CALL SLACPY( 'A', M, 1, A, LDA, U, LDU )
00559 *           computing all M left singular vectors of the M x 1 matrix
00560             IF ( N1 .NE. N  ) THEN
00561                CALL SGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR )
00562                CALL SORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR )
00563                CALL SCOPY( M, A(1,1), 1, U(1,1), 1 )
00564             END IF
00565          END IF
00566          IF ( RSVEC ) THEN
00567              V(1,1) = ONE
00568          END IF
00569          IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
00570             SVA(1)  = SVA(1) / SCALEM
00571             SCALEM  = ONE
00572          END IF
00573          WORK(1) = ONE / SCALEM
00574          WORK(2) = ONE
00575          IF ( SVA(1) .NE. ZERO ) THEN
00576             IWORK(1) = 1
00577             IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
00578                IWORK(2) = 1
00579             ELSE
00580                IWORK(2) = 0
00581             END IF
00582          ELSE
00583             IWORK(1) = 0
00584             IWORK(2) = 0
00585          END IF
00586          IF ( ERREST ) WORK(3) = ONE
00587          IF ( LSVEC .AND. RSVEC ) THEN
00588             WORK(4) = ONE
00589             WORK(5) = ONE
00590          END IF
00591          IF ( L2TRAN ) THEN
00592             WORK(6) = ZERO
00593             WORK(7) = ZERO
00594          END IF
00595          RETURN
00596 *
00597       END IF
00598 *
00599       TRANSP = .FALSE.
00600       L2TRAN = L2TRAN .AND. ( M .EQ. N )
00601 *
00602       AATMAX = -ONE
00603       AATMIN =  BIG
00604       IF ( ROWPIV .OR. L2TRAN ) THEN
00605 *
00606 *     Compute the row norms, needed to determine row pivoting sequence
00607 *     (in the case of heavily row weighted A, row pivoting is strongly
00608 *     advised) and to collect information needed to compare the
00609 *     structures of A * A^t and A^t * A (in the case L2TRAN.EQ..TRUE.).
00610 *
00611          IF ( L2TRAN ) THEN
00612             DO 1950 p = 1, M
00613                XSC   = ZERO
00614                TEMP1 = ONE
00615                CALL SLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
00616 *              SLASSQ gets both the ell_2 and the ell_infinity norm
00617 *              in one pass through the vector
00618                WORK(M+N+p)  = XSC * SCALEM
00619                WORK(N+p)    = XSC * (SCALEM*SQRT(TEMP1))
00620                AATMAX = AMAX1( AATMAX, WORK(N+p) )
00621                IF (WORK(N+p) .NE. ZERO) AATMIN = AMIN1(AATMIN,WORK(N+p))
00622  1950       CONTINUE
00623          ELSE
00624             DO 1904 p = 1, M
00625                WORK(M+N+p) = SCALEM*ABS( A(p,ISAMAX(N,A(p,1),LDA)) )
00626                AATMAX = AMAX1( AATMAX, WORK(M+N+p) )
00627                AATMIN = AMIN1( AATMIN, WORK(M+N+p) )
00628  1904       CONTINUE
00629          END IF
00630 *
00631       END IF
00632 *
00633 *     For square matrix A try to determine whether A^t  would be  better
00634 *     input for the preconditioned Jacobi SVD, with faster convergence.
00635 *     The decision is based on an O(N) function of the vector of column
00636 *     and row norms of A, based on the Shannon entropy. This should give
00637 *     the right choice in most cases when the difference actually matters.
00638 *     It may fail and pick the slower converging side.
00639 *
00640       ENTRA  = ZERO
00641       ENTRAT = ZERO
00642       IF ( L2TRAN ) THEN
00643 *
00644          XSC   = ZERO
00645          TEMP1 = ONE
00646          CALL SLASSQ( N, SVA, 1, XSC, TEMP1 )
00647          TEMP1 = ONE / TEMP1
00648 *
00649          ENTRA = ZERO
00650          DO 1113 p = 1, N
00651             BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1
00652             IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * ALOG(BIG1)
00653  1113    CONTINUE
00654          ENTRA = - ENTRA / ALOG(FLOAT(N))
00655 *
00656 *        Now, SVA().^2/Trace(A^t * A) is a point in the probability simplex.
00657 *        It is derived from the diagonal of  A^t * A.  Do the same with the
00658 *        diagonal of A * A^t, compute the entropy of the corresponding
00659 *        probability distribution. Note that A * A^t and A^t * A have the
00660 *        same trace.
00661 *
00662          ENTRAT = ZERO
00663          DO 1114 p = N+1, N+M
00664             BIG1 = ( ( WORK(p) / XSC )**2 ) * TEMP1
00665             IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * ALOG(BIG1)
00666  1114    CONTINUE
00667          ENTRAT = - ENTRAT / ALOG(FLOAT(M))
00668 *
00669 *        Analyze the entropies and decide A or A^t. Smaller entropy
00670 *        usually means better input for the algorithm.
00671 *
00672          TRANSP = ( ENTRAT .LT. ENTRA )
00673 *
00674 *        If A^t is better than A, transpose A.
00675 *
00676          IF ( TRANSP ) THEN
00677 *           In an optimal implementation, this trivial transpose
00678 *           should be replaced with faster transpose.
00679             DO 1115 p = 1, N - 1
00680                DO 1116 q = p + 1, N
00681                    TEMP1 = A(q,p)
00682                   A(q,p) = A(p,q)
00683                   A(p,q) = TEMP1
00684  1116          CONTINUE
00685  1115       CONTINUE
00686             DO 1117 p = 1, N
00687                WORK(M+N+p) = SVA(p)
00688                SVA(p)      = WORK(N+p)
00689  1117       CONTINUE
00690             TEMP1  = AAPP
00691             AAPP   = AATMAX
00692             AATMAX = TEMP1
00693             TEMP1  = AAQQ
00694             AAQQ   = AATMIN
00695             AATMIN = TEMP1
00696             KILL   = LSVEC
00697             LSVEC  = RSVEC
00698             RSVEC  = KILL
00699             IF ( LSVEC ) N1 = N 
00700 *
00701             ROWPIV = .TRUE.
00702          END IF
00703 *
00704       END IF
00705 *     END IF L2TRAN
00706 *
00707 *     Scale the matrix so that its maximal singular value remains less
00708 *     than SQRT(BIG) -- the matrix is scaled so that its maximal column
00709 *     has Euclidean norm equal to SQRT(BIG/N). The only reason to keep
00710 *     SQRT(BIG) instead of BIG is the fact that SGEJSV uses LAPACK and
00711 *     BLAS routines that, in some implementations, are not capable of
00712 *     working in the full interval [SFMIN,BIG] and that they may provoke
00713 *     overflows in the intermediate results. If the singular values spread
00714 *     from SFMIN to BIG, then SGESVJ will compute them. So, in that case,
00715 *     one should use SGESVJ instead of SGEJSV.
00716 *
00717       BIG1   = SQRT( BIG )
00718       TEMP1  = SQRT( BIG / FLOAT(N) )
00719 *
00720       CALL SLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
00721       IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
00722           AAQQ = ( AAQQ / AAPP ) * TEMP1
00723       ELSE
00724           AAQQ = ( AAQQ * TEMP1 ) / AAPP
00725       END IF
00726       TEMP1 = TEMP1 * SCALEM
00727       CALL SLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
00728 *
00729 *     To undo scaling at the end of this procedure, multiply the
00730 *     computed singular values with USCAL2 / USCAL1.
00731 *
00732       USCAL1 = TEMP1
00733       USCAL2 = AAPP
00734 *
00735       IF ( L2KILL ) THEN
00736 *        L2KILL enforces computation of nonzero singular values in
00737 *        the restricted range of condition number of the initial A,
00738 *        sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN).
00739          XSC = SQRT( SFMIN )
00740       ELSE
00741          XSC = SMALL
00742 *
00743 *        Now, if the condition number of A is too big,
00744 *        sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN,
00745 *        as a precaution measure, the full SVD is computed using SGESVJ
00746 *        with accumulated Jacobi rotations. This provides numerically
00747 *        more robust computation, at the cost of slightly increased run
00748 *        time. Depending on the concrete implementation of BLAS and LAPACK
00749 *        (i.e. how they behave in presence of extreme ill-conditioning) the
00750 *        implementor may decide to remove this switch.
00751          IF ( ( AAQQ.LT.SQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
00752             JRACC = .TRUE.
00753          END IF
00754 *
00755       END IF
00756       IF ( AAQQ .LT. XSC ) THEN
00757          DO 700 p = 1, N
00758             IF ( SVA(p) .LT. XSC ) THEN
00759                CALL SLASET( 'A', M, 1, ZERO, ZERO, A(1,p), LDA )
00760                SVA(p) = ZERO
00761             END IF
00762  700     CONTINUE
00763       END IF
00764 *
00765 *     Preconditioning using QR factorization with pivoting
00766 *
00767       IF ( ROWPIV ) THEN
00768 *        Optional row permutation (Bjoerck row pivoting):
00769 *        A result by Cox and Higham shows that the Bjoerck's
00770 *        row pivoting combined with standard column pivoting
00771 *        has similar effect as Powell-Reid complete pivoting.
00772 *        The ell-infinity norms of A are made nonincreasing.
00773          DO 1952 p = 1, M - 1
00774             q = ISAMAX( M-p+1, WORK(M+N+p), 1 ) + p - 1
00775             IWORK(2*N+p) = q
00776             IF ( p .NE. q ) THEN
00777                TEMP1       = WORK(M+N+p)
00778                WORK(M+N+p) = WORK(M+N+q)
00779                WORK(M+N+q) = TEMP1
00780             END IF
00781  1952    CONTINUE
00782          CALL SLASWP( N, A, LDA, 1, M-1, IWORK(2*N+1), 1 )
00783       END IF
00784 *
00785 *     End of the preparation phase (scaling, optional sorting and
00786 *     transposing, optional flushing of small columns).
00787 *
00788 *     Preconditioning
00789 *
00790 *     If the full SVD is needed, the right singular vectors are computed
00791 *     from a matrix equation, and for that we need theoretical analysis
00792 *     of the Businger-Golub pivoting. So we use SGEQP3 as the first RR QRF.
00793 *     In all other cases the first RR QRF can be chosen by other criteria
00794 *     (eg speed by replacing global with restricted window pivoting, such
00795 *     as in SGEQPX from TOMS # 782). Good results will be obtained using
00796 *     SGEQPX with properly (!) chosen numerical parameters.
00797 *     Any improvement of SGEQP3 improves overal performance of SGEJSV.
00798 *
00799 *     A * P1 = Q1 * [ R1^t 0]^t:
00800       DO 1963 p = 1, N
00801 *        .. all columns are free columns
00802          IWORK(p) = 0
00803  1963 CONTINUE
00804       CALL SGEQP3( M,N,A,LDA, IWORK,WORK, WORK(N+1),LWORK-N, IERR )
00805 *
00806 *     The upper triangular matrix R1 from the first QRF is inspected for
00807 *     rank deficiency and possibilities for deflation, or possible
00808 *     ill-conditioning. Depending on the user specified flag L2RANK,
00809 *     the procedure explores possibilities to reduce the numerical
00810 *     rank by inspecting the computed upper triangular factor. If
00811 *     L2RANK or L2ABER are up, then SGEJSV will compute the SVD of
00812 *     A + dA, where ||dA|| <= f(M,N)*EPSLN.
00813 *
00814       NR = 1
00815       IF ( L2ABER ) THEN
00816 *        Standard absolute error bound suffices. All sigma_i with
00817 *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
00818 *        agressive enforcement of lower numerical rank by introducing a
00819 *        backward error of the order of N*EPSLN*||A||.
00820          TEMP1 = SQRT(FLOAT(N))*EPSLN
00821          DO 3001 p = 2, N
00822             IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN
00823                NR = NR + 1
00824             ELSE
00825                GO TO 3002
00826             END IF
00827  3001    CONTINUE
00828  3002    CONTINUE
00829       ELSE IF ( L2RANK ) THEN
00830 *        .. similarly as above, only slightly more gentle (less agressive).
00831 *        Sudden drop on the diagonal of R1 is used as the criterion for
00832 *        close-to-rank-defficient.
00833          TEMP1 = SQRT(SFMIN)
00834          DO 3401 p = 2, N
00835             IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
00836      &           ( ABS(A(p,p)) .LT. SMALL ) .OR.
00837      &           ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
00838             NR = NR + 1
00839  3401    CONTINUE
00840  3402    CONTINUE
00841 *
00842       ELSE
00843 *        The goal is high relative accuracy. However, if the matrix
00844 *        has high scaled condition number the relative accuracy is in
00845 *        general not feasible. Later on, a condition number estimator
00846 *        will be deployed to estimate the scaled condition number.
00847 *        Here we just remove the underflowed part of the triangular
00848 *        factor. This prevents the situation in which the code is
00849 *        working hard to get the accuracy not warranted by the data.
00850          TEMP1  = SQRT(SFMIN)
00851          DO 3301 p = 2, N
00852             IF ( ( ABS(A(p,p)) .LT. SMALL ) .OR.
00853      &           ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
00854             NR = NR + 1
00855  3301    CONTINUE
00856  3302    CONTINUE
00857 *
00858       END IF
00859 *
00860       ALMORT = .FALSE.
00861       IF ( NR .EQ. N ) THEN
00862          MAXPRJ = ONE
00863          DO 3051 p = 2, N
00864             TEMP1  = ABS(A(p,p)) / SVA(IWORK(p))
00865             MAXPRJ = AMIN1( MAXPRJ, TEMP1 )
00866  3051    CONTINUE
00867          IF ( MAXPRJ**2 .GE. ONE - FLOAT(N)*EPSLN ) ALMORT = .TRUE.
00868       END IF
00869 *
00870 *
00871       SCONDA = - ONE
00872       CONDR1 = - ONE
00873       CONDR2 = - ONE
00874 *
00875       IF ( ERREST ) THEN
00876          IF ( N .EQ. NR ) THEN
00877             IF ( RSVEC ) THEN
00878 *              .. V is available as workspace
00879                CALL SLACPY( 'U', N, N, A, LDA, V, LDV )
00880                DO 3053 p = 1, N
00881                   TEMP1 = SVA(IWORK(p))
00882                   CALL SSCAL( p, ONE/TEMP1, V(1,p), 1 )
00883  3053          CONTINUE
00884                CALL SPOCON( 'U', N, V, LDV, ONE, TEMP1,
00885      &              WORK(N+1), IWORK(2*N+M+1), IERR )
00886             ELSE IF ( LSVEC ) THEN
00887 *              .. U is available as workspace
00888                CALL SLACPY( 'U', N, N, A, LDA, U, LDU )
00889                DO 3054 p = 1, N
00890                   TEMP1 = SVA(IWORK(p))
00891                   CALL SSCAL( p, ONE/TEMP1, U(1,p), 1 )
00892  3054          CONTINUE
00893                CALL SPOCON( 'U', N, U, LDU, ONE, TEMP1,
00894      &              WORK(N+1), IWORK(2*N+M+1), IERR )
00895             ELSE
00896                CALL SLACPY( 'U', N, N, A, LDA, WORK(N+1), N )
00897                DO 3052 p = 1, N
00898                   TEMP1 = SVA(IWORK(p))
00899                   CALL SSCAL( p, ONE/TEMP1, WORK(N+(p-1)*N+1), 1 )
00900  3052          CONTINUE
00901 *           .. the columns of R are scaled to have unit Euclidean lengths.
00902                CALL SPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1,
00903      &              WORK(N+N*N+1), IWORK(2*N+M+1), IERR )
00904             END IF
00905             SCONDA = ONE / SQRT(TEMP1)
00906 *           SCONDA is an estimate of SQRT(||(R^t * R)^(-1)||_1).
00907 *           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
00908          ELSE
00909             SCONDA = - ONE
00910          END IF
00911       END IF
00912 *
00913       L2PERT = L2PERT .AND. ( ABS( A(1,1)/A(NR,NR) ) .GT. SQRT(BIG1) )
00914 *     If there is no violent scaling, artificial perturbation is not needed.
00915 *
00916 *     Phase 3:
00917 *
00918       IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
00919 *
00920 *         Singular Values only
00921 *
00922 *         .. transpose A(1:NR,1:N)
00923          DO 1946 p = 1, MIN0( N-1, NR )
00924             CALL SCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
00925  1946    CONTINUE
00926 *
00927 *        The following two DO-loops introduce small relative perturbation
00928 *        into the strict upper triangle of the lower triangular matrix.
00929 *        Small entries below the main diagonal are also changed.
00930 *        This modification is useful if the computing environment does not
00931 *        provide/allow FLUSH TO ZERO underflow, for it prevents many
00932 *        annoying denormalized numbers in case of strongly scaled matrices.
00933 *        The perturbation is structured so that it does not introduce any
00934 *        new perturbation of the singular values, and it does not destroy
00935 *        the job done by the preconditioner.
00936 *        The licence for this perturbation is in the variable L2PERT, which
00937 *        should be .FALSE. if FLUSH TO ZERO underflow is active.
00938 *
00939          IF ( .NOT. ALMORT ) THEN
00940 *
00941             IF ( L2PERT ) THEN
00942 *              XSC = SQRT(SMALL)
00943                XSC = EPSLN / FLOAT(N)
00944                DO 4947 q = 1, NR
00945                   TEMP1 = XSC*ABS(A(q,q))
00946                   DO 4949 p = 1, N
00947                      IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
00948      &                    .OR. ( p .LT. q ) )
00949      &                     A(p,q) = SIGN( TEMP1, A(p,q) )
00950  4949             CONTINUE
00951  4947          CONTINUE
00952             ELSE
00953                CALL SLASET( 'U', NR-1,NR-1, ZERO,ZERO, A(1,2),LDA )
00954             END IF
00955 *
00956 *            .. second preconditioning using the QR factorization
00957 *
00958             CALL SGEQRF( N,NR, A,LDA, WORK, WORK(N+1),LWORK-N, IERR )
00959 *
00960 *           .. and transpose upper to lower triangular
00961             DO 1948 p = 1, NR - 1
00962                CALL SCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
00963  1948       CONTINUE
00964 *
00965          END IF
00966 *
00967 *           Row-cyclic Jacobi SVD algorithm with column pivoting
00968 *
00969 *           .. again some perturbation (a "background noise") is added
00970 *           to drown denormals
00971             IF ( L2PERT ) THEN
00972 *              XSC = SQRT(SMALL)
00973                XSC = EPSLN / FLOAT(N)
00974                DO 1947 q = 1, NR
00975                   TEMP1 = XSC*ABS(A(q,q))
00976                   DO 1949 p = 1, NR
00977                      IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
00978      &                       .OR. ( p .LT. q ) )
00979      &                   A(p,q) = SIGN( TEMP1, A(p,q) )
00980  1949             CONTINUE
00981  1947          CONTINUE
00982             ELSE
00983                CALL SLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA )
00984             END IF
00985 *
00986 *           .. and one-sided Jacobi rotations are started on a lower
00987 *           triangular matrix (plus perturbation which is ignored in
00988 *           the part which destroys triangular form (confusing?!))
00989 *
00990             CALL SGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA,
00991      &                      N, V, LDV, WORK, LWORK, INFO )
00992 *
00993             SCALEM  = WORK(1)
00994             NUMRANK = NINT(WORK(2))
00995 *
00996 *
00997       ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
00998 *
00999 *        -> Singular Values and Right Singular Vectors <-
01000 *
01001          IF ( ALMORT ) THEN
01002 *
01003 *           .. in this case NR equals N
01004             DO 1998 p = 1, NR
01005                CALL SCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
01006  1998       CONTINUE
01007             CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01008 *
01009             CALL SGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA,
01010      &                  WORK, LWORK, INFO )
01011             SCALEM  = WORK(1)
01012             NUMRANK = NINT(WORK(2))
01013 
01014          ELSE
01015 *
01016 *        .. two more QR factorizations ( one QRF is not enough, two require
01017 *        accumulated product of Jacobi rotations, three are perfect )
01018 *
01019             CALL SLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA )
01020             CALL SGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR)
01021             CALL SLACPY( 'Lower', NR, NR, A, LDA, V, LDV )
01022             CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01023             CALL SGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1),
01024      &                   LWORK-2*N, IERR )
01025             DO 8998 p = 1, NR
01026                CALL SCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
01027  8998       CONTINUE
01028             CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01029 *
01030             CALL SGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
01031      &                  LDU, WORK(N+1), LWORK, INFO )
01032             SCALEM  = WORK(N+1)
01033             NUMRANK = NINT(WORK(N+2))
01034             IF ( NR .LT. N ) THEN
01035                CALL SLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1),   LDV )
01036                CALL SLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1),   LDV )
01037                CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV )
01038             END IF
01039 *
01040          CALL SORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK,
01041      &               V, LDV, WORK(N+1), LWORK-N, IERR )
01042 *
01043          END IF
01044 *
01045          DO 8991 p = 1, N
01046             CALL SCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
01047  8991    CONTINUE
01048          CALL SLACPY( 'All', N, N, A, LDA, V, LDV )
01049 *
01050          IF ( TRANSP ) THEN
01051             CALL SLACPY( 'All', N, N, V, LDV, U, LDU )
01052          END IF
01053 *
01054       ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
01055 *
01056 *        .. Singular Values and Left Singular Vectors                 ..
01057 *
01058 *        .. second preconditioning step to avoid need to accumulate
01059 *        Jacobi rotations in the Jacobi iterations.
01060          DO 1965 p = 1, NR
01061             CALL SCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
01062  1965    CONTINUE
01063          CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
01064 *
01065          CALL SGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1),
01066      &              LWORK-2*N, IERR )
01067 *
01068          DO 1967 p = 1, NR - 1
01069             CALL SCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
01070  1967    CONTINUE
01071          CALL SLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
01072 *
01073          CALL SGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
01074      &        LDA, WORK(N+1), LWORK-N, INFO )
01075          SCALEM  = WORK(N+1)
01076          NUMRANK = NINT(WORK(N+2))
01077 *
01078          IF ( NR .LT. M ) THEN
01079             CALL SLASET( 'A',  M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU )
01080             IF ( NR .LT. N1 ) THEN
01081                CALL SLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU )
01082                CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU )
01083             END IF
01084          END IF
01085 *
01086          CALL SORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
01087      &               LDU, WORK(N+1), LWORK-N, IERR )
01088 *
01089          IF ( ROWPIV )
01090      &       CALL SLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
01091 *
01092          DO 1974 p = 1, N1
01093             XSC = ONE / SNRM2( M, U(1,p), 1 )
01094             CALL SSCAL( M, XSC, U(1,p), 1 )
01095  1974    CONTINUE
01096 *
01097          IF ( TRANSP ) THEN
01098             CALL SLACPY( 'All', N, N, U, LDU, V, LDV )
01099          END IF
01100 *
01101       ELSE
01102 *
01103 *        .. Full SVD ..
01104 *
01105          IF ( .NOT. JRACC ) THEN
01106 *
01107          IF ( .NOT. ALMORT ) THEN
01108 *
01109 *           Second Preconditioning Step (QRF [with pivoting])
01110 *           Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
01111 *           equivalent to an LQF CALL. Since in many libraries the QRF
01112 *           seems to be better optimized than the LQF, we do explicit
01113 *           transpose and use the QRF. This is subject to changes in an
01114 *           optimized implementation of SGEJSV.
01115 *
01116             DO 1968 p = 1, NR
01117                CALL SCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
01118  1968       CONTINUE
01119 *
01120 *           .. the following two loops perturb small entries to avoid
01121 *           denormals in the second QR factorization, where they are
01122 *           as good as zeros. This is done to avoid painfully slow
01123 *           computation with denormals. The relative size of the perturbation
01124 *           is a parameter that can be changed by the implementer.
01125 *           This perturbation device will be obsolete on machines with
01126 *           properly implemented arithmetic.
01127 *           To switch it off, set L2PERT=.FALSE. To remove it from  the
01128 *           code, remove the action under L2PERT=.TRUE., leave the ELSE part.
01129 *           The following two loops should be blocked and fused with the
01130 *           transposed copy above.
01131 *
01132             IF ( L2PERT ) THEN
01133                XSC = SQRT(SMALL)
01134                DO 2969 q = 1, NR
01135                   TEMP1 = XSC*ABS( V(q,q) )
01136                   DO 2968 p = 1, N
01137                      IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
01138      &                   .OR. ( p .LT. q ) )
01139      &                   V(p,q) = SIGN( TEMP1, V(p,q) )
01140                      IF ( p. LT. q ) V(p,q) = - V(p,q)
01141  2968             CONTINUE
01142  2969          CONTINUE
01143             ELSE
01144                CALL SLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01145             END IF
01146 *
01147 *           Estimate the row scaled condition number of R1
01148 *           (If R1 is rectangular, N > NR, then the condition number
01149 *           of the leading NR x NR submatrix is estimated.)
01150 *
01151             CALL SLACPY( 'L', NR, NR, V, LDV, WORK(2*N+1), NR )
01152             DO 3950 p = 1, NR
01153                TEMP1 = SNRM2(NR-p+1,WORK(2*N+(p-1)*NR+p),1)
01154                CALL SSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1)
01155  3950       CONTINUE
01156             CALL SPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
01157      &                   WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
01158             CONDR1 = ONE / SQRT(TEMP1)
01159 *           .. here need a second oppinion on the condition number
01160 *           .. then assume worst case scenario
01161 *           R1 is OK for inverse <=> CONDR1 .LT. FLOAT(N)
01162 *           more conservative    <=> CONDR1 .LT. SQRT(FLOAT(N))
01163 *
01164             COND_OK = SQRT(FLOAT(NR))
01165 *[TP]       COND_OK is a tuning parameter.
01166 
01167             IF ( CONDR1 .LT. COND_OK ) THEN
01168 *              .. the second QRF without pivoting. Note: in an optimized
01169 *              implementation, this QRF should be implemented as the QRF
01170 *              of a lower triangular matrix.
01171 *              R1^t = Q2 * R2
01172                CALL SGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
01173      &              LWORK-2*N, IERR )
01174 *
01175                IF ( L2PERT ) THEN
01176                   XSC = SQRT(SMALL)/EPSLN
01177                   DO 3959 p = 2, NR
01178                      DO 3958 q = 1, p - 1
01179                         TEMP1 = XSC * AMIN1(ABS(V(p,p)),ABS(V(q,q)))
01180                         IF ( ABS(V(q,p)) .LE. TEMP1 )
01181      &                     V(q,p) = SIGN( TEMP1, V(q,p) )
01182  3958                CONTINUE
01183  3959             CONTINUE
01184                END IF
01185 *
01186                IF ( NR .NE. N )
01187 *              .. save ...
01188      &         CALL SLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
01189 *
01190 *           .. this transposed copy should be better than naive
01191                DO 1969 p = 1, NR - 1
01192                   CALL SCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
01193  1969          CONTINUE
01194 *
01195                CONDR2 = CONDR1
01196 *
01197             ELSE
01198 *
01199 *              .. ill-conditioned case: second QRF with pivoting
01200 *              Note that windowed pivoting would be equaly good
01201 *              numerically, and more run-time efficient. So, in
01202 *              an optimal implementation, the next call to SGEQP3
01203 *              should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
01204 *              with properly (carefully) chosen parameters.
01205 *
01206 *              R1^t * P2 = Q2 * R2
01207                DO 3003 p = 1, NR
01208                   IWORK(N+p) = 0
01209  3003          CONTINUE
01210                CALL SGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1),
01211      &                  WORK(2*N+1), LWORK-2*N, IERR )
01212 **               CALL SGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
01213 **     &              LWORK-2*N, IERR )
01214                IF ( L2PERT ) THEN
01215                   XSC = SQRT(SMALL)
01216                   DO 3969 p = 2, NR
01217                      DO 3968 q = 1, p - 1
01218                         TEMP1 = XSC * AMIN1(ABS(V(p,p)),ABS(V(q,q)))
01219                         IF ( ABS(V(q,p)) .LE. TEMP1 )
01220      &                     V(q,p) = SIGN( TEMP1, V(q,p) )
01221  3968                CONTINUE
01222  3969             CONTINUE
01223                END IF
01224 *
01225                CALL SLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
01226 *
01227                IF ( L2PERT ) THEN
01228                   XSC = SQRT(SMALL)
01229                   DO 8970 p = 2, NR
01230                      DO 8971 q = 1, p - 1
01231                         TEMP1 = XSC * AMIN1(ABS(V(p,p)),ABS(V(q,q)))
01232                         V(p,q) = - SIGN( TEMP1, V(q,p) )
01233  8971                CONTINUE
01234  8970             CONTINUE
01235                ELSE
01236                   CALL SLASET( 'L',NR-1,NR-1,ZERO,ZERO,V(2,1),LDV )
01237                END IF
01238 *              Now, compute R2 = L3 * Q3, the LQ factorization.
01239                CALL SGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1),
01240      &               WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
01241 *              .. and estimate the condition number
01242                CALL SLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR )
01243                DO 4950 p = 1, NR
01244                   TEMP1 = SNRM2( p, WORK(2*N+N*NR+NR+p), NR )
01245                   CALL SSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR )
01246  4950          CONTINUE
01247                CALL SPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
01248      &              WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR )
01249                CONDR2 = ONE / SQRT(TEMP1)
01250 *
01251                IF ( CONDR2 .GE. COND_OK ) THEN
01252 *                 .. save the Householder vectors used for Q3
01253 *                 (this overwrittes the copy of R2, as it will not be
01254 *                 needed in this branch, but it does not overwritte the
01255 *                 Huseholder vectors of Q2.).
01256                   CALL SLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
01257 *                 .. and the rest of the information on Q3 is in
01258 *                 WORK(2*N+N*NR+1:2*N+N*NR+N)
01259                END IF
01260 *
01261             END IF
01262 *
01263             IF ( L2PERT ) THEN
01264                XSC = SQRT(SMALL)
01265                DO 4968 q = 2, NR
01266                   TEMP1 = XSC * V(q,q)
01267                   DO 4969 p = 1, q - 1
01268 *                    V(p,q) = - SIGN( TEMP1, V(q,p) )
01269                      V(p,q) = - SIGN( TEMP1, V(p,q) )
01270  4969             CONTINUE
01271  4968          CONTINUE
01272             ELSE
01273                CALL SLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
01274             END IF
01275 *
01276 *        Second preconditioning finished; continue with Jacobi SVD
01277 *        The input matrix is lower trinagular.
01278 *
01279 *        Recover the right singular vectors as solution of a well
01280 *        conditioned triangular matrix equation.
01281 *
01282             IF ( CONDR1 .LT. COND_OK ) THEN
01283 *
01284                CALL SGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U,
01285      &              LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO )
01286                SCALEM  = WORK(2*N+N*NR+NR+1)
01287                NUMRANK = NINT(WORK(2*N+N*NR+NR+2))
01288                DO 3970 p = 1, NR
01289                   CALL SCOPY( NR, V(1,p), 1, U(1,p), 1 )
01290                   CALL SSCAL( NR, SVA(p),    V(1,p), 1 )
01291  3970          CONTINUE
01292 
01293 *        .. pick the right matrix equation and solve it
01294 *
01295                IF ( NR. EQ. N ) THEN
01296 * :))             .. best case, R1 is inverted. The solution of this matrix
01297 *                 equation is Q2*V2 = the product of the Jacobi rotations
01298 *                 used in SGESVJ, premultiplied with the orthogonal matrix
01299 *                 from the second QR factorization.
01300                   CALL STRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV )
01301                ELSE
01302 *                 .. R1 is well conditioned, but non-square. Transpose(R2)
01303 *                 is inverted to get the product of the Jacobi rotations
01304 *                 used in SGESVJ. The Q-factor from the second QR
01305 *                 factorization is then built in explicitly.
01306                   CALL STRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1),
01307      &                 N,V,LDV)
01308                   IF ( NR .LT. N ) THEN
01309                     CALL SLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
01310                     CALL SLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
01311                     CALL SLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
01312                   END IF
01313                   CALL SORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
01314      &                 V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
01315                END IF
01316 *
01317             ELSE IF ( CONDR2 .LT. COND_OK ) THEN
01318 *
01319 * :)           .. the input matrix A is very likely a relative of
01320 *              the Kahan matrix :)
01321 *              The matrix R2 is inverted. The solution of the matrix equation
01322 *              is Q3^T*V3 = the product of the Jacobi rotations (appplied to
01323 *              the lower triangular L3 from the LQ factorization of
01324 *              R2=L3*Q3), pre-multiplied with the transposed Q3.
01325                CALL SGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
01326      &              LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
01327                SCALEM  = WORK(2*N+N*NR+NR+1)
01328                NUMRANK = NINT(WORK(2*N+N*NR+NR+2))
01329                DO 3870 p = 1, NR
01330                   CALL SCOPY( NR, V(1,p), 1, U(1,p), 1 )
01331                   CALL SSCAL( NR, SVA(p),    U(1,p), 1 )
01332  3870          CONTINUE
01333                CALL STRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU)
01334 *              .. apply the permutation from the second QR factorization
01335                DO 873 q = 1, NR
01336                   DO 872 p = 1, NR
01337                      WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
01338  872              CONTINUE
01339                   DO 874 p = 1, NR
01340                      U(p,q) = WORK(2*N+N*NR+NR+p)
01341  874              CONTINUE
01342  873           CONTINUE
01343                IF ( NR .LT. N ) THEN
01344                   CALL SLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
01345                   CALL SLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
01346                   CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
01347                END IF
01348                CALL SORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
01349      &              V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
01350             ELSE
01351 *              Last line of defense.
01352 * #:(          This is a rather pathological case: no scaled condition
01353 *              improvement after two pivoted QR factorizations. Other
01354 *              possibility is that the rank revealing QR factorization
01355 *              or the condition estimator has failed, or the COND_OK
01356 *              is set very close to ONE (which is unnecessary). Normally,
01357 *              this branch should never be executed, but in rare cases of
01358 *              failure of the RRQR or condition estimator, the last line of
01359 *              defense ensures that SGEJSV completes the task.
01360 *              Compute the full SVD of L3 using SGESVJ with explicit
01361 *              accumulation of Jacobi rotations.
01362                CALL SGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
01363      &              LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
01364                SCALEM  = WORK(2*N+N*NR+NR+1)
01365                NUMRANK = NINT(WORK(2*N+N*NR+NR+2))
01366                IF ( NR .LT. N ) THEN
01367                   CALL SLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
01368                   CALL SLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
01369                   CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
01370                END IF
01371                CALL SORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
01372      &              V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
01373 *
01374                CALL SORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N,
01375      &              WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1),
01376      &              LWORK-2*N-N*NR-NR, IERR )
01377                DO 773 q = 1, NR
01378                   DO 772 p = 1, NR
01379                      WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
01380  772              CONTINUE
01381                   DO 774 p = 1, NR
01382                      U(p,q) = WORK(2*N+N*NR+NR+p)
01383  774              CONTINUE
01384  773           CONTINUE
01385 *
01386             END IF
01387 *
01388 *           Permute the rows of V using the (column) permutation from the
01389 *           first QRF. Also, scale the columns to make them unit in
01390 *           Euclidean norm. This applies to all cases.
01391 *
01392             TEMP1 = SQRT(FLOAT(N)) * EPSLN
01393             DO 1972 q = 1, N
01394                DO 972 p = 1, N
01395                   WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
01396   972          CONTINUE
01397                DO 973 p = 1, N
01398                   V(p,q) = WORK(2*N+N*NR+NR+p)
01399   973          CONTINUE
01400                XSC = ONE / SNRM2( N, V(1,q), 1 )
01401                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01402      &           CALL SSCAL( N, XSC, V(1,q), 1 )
01403  1972       CONTINUE
01404 *           At this moment, V contains the right singular vectors of A.
01405 *           Next, assemble the left singular vector matrix U (M x N).
01406             IF ( NR .LT. M ) THEN
01407                CALL SLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
01408                IF ( NR .LT. N1 ) THEN
01409                   CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
01410                   CALL SLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU)
01411                END IF
01412             END IF
01413 *
01414 *           The Q matrix from the first QRF is built into the left singular
01415 *           matrix U. This applies to all cases.
01416 *
01417             CALL SORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U,
01418      &           LDU, WORK(N+1), LWORK-N, IERR )
01419 
01420 *           The columns of U are normalized. The cost is O(M*N) flops.
01421             TEMP1 = SQRT(FLOAT(M)) * EPSLN
01422             DO 1973 p = 1, NR
01423                XSC = ONE / SNRM2( M, U(1,p), 1 )
01424                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01425      &          CALL SSCAL( M, XSC, U(1,p), 1 )
01426  1973       CONTINUE
01427 *
01428 *           If the initial QRF is computed with row pivoting, the left
01429 *           singular vectors must be adjusted.
01430 *
01431             IF ( ROWPIV )
01432      &          CALL SLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
01433 *
01434          ELSE
01435 *
01436 *        .. the initial matrix A has almost orthogonal columns and
01437 *        the second QRF is not needed
01438 *
01439             CALL SLACPY( 'Upper', N, N, A, LDA, WORK(N+1), N )
01440             IF ( L2PERT ) THEN
01441                XSC = SQRT(SMALL)
01442                DO 5970 p = 2, N
01443                   TEMP1 = XSC * WORK( N + (p-1)*N + p )
01444                   DO 5971 q = 1, p - 1
01445                      WORK(N+(q-1)*N+p)=-SIGN(TEMP1,WORK(N+(p-1)*N+q))
01446  5971             CONTINUE
01447  5970          CONTINUE
01448             ELSE
01449                CALL SLASET( 'Lower',N-1,N-1,ZERO,ZERO,WORK(N+2),N )
01450             END IF
01451 *
01452             CALL SGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA,
01453      &           N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO )
01454 *
01455             SCALEM  = WORK(N+N*N+1)
01456             NUMRANK = NINT(WORK(N+N*N+2))
01457             DO 6970 p = 1, N
01458                CALL SCOPY( N, WORK(N+(p-1)*N+1), 1, U(1,p), 1 )
01459                CALL SSCAL( N, SVA(p), WORK(N+(p-1)*N+1), 1 )
01460  6970       CONTINUE
01461 *
01462             CALL STRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N,
01463      &           ONE, A, LDA, WORK(N+1), N )
01464             DO 6972 p = 1, N
01465                CALL SCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV )
01466  6972       CONTINUE
01467             TEMP1 = SQRT(FLOAT(N))*EPSLN
01468             DO 6971 p = 1, N
01469                XSC = ONE / SNRM2( N, V(1,p), 1 )
01470                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01471      &            CALL SSCAL( N, XSC, V(1,p), 1 )
01472  6971       CONTINUE
01473 *
01474 *           Assemble the left singular vector matrix U (M x N).
01475 *
01476             IF ( N .LT. M ) THEN
01477                CALL SLASET( 'A',  M-N, N, ZERO, ZERO, U(N+1,1), LDU )
01478                IF ( N .LT. N1 ) THEN
01479                   CALL SLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
01480                   CALL SLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
01481                END IF
01482             END IF
01483             CALL SORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
01484      &           LDU, WORK(N+1), LWORK-N, IERR )
01485             TEMP1 = SQRT(FLOAT(M))*EPSLN
01486             DO 6973 p = 1, N1
01487                XSC = ONE / SNRM2( M, U(1,p), 1 )
01488                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01489      &            CALL SSCAL( M, XSC, U(1,p), 1 )
01490  6973       CONTINUE
01491 *
01492             IF ( ROWPIV )
01493      &         CALL SLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
01494 *
01495          END IF
01496 *
01497 *        end of the  >> almost orthogonal case <<  in the full SVD
01498 *
01499          ELSE
01500 *
01501 *        This branch deploys a preconditioned Jacobi SVD with explicitly
01502 *        accumulated rotations. It is included as optional, mainly for
01503 *        experimental purposes. It does perfom well, and can also be used.
01504 *        In this implementation, this branch will be automatically activated
01505 *        if the  condition number sigma_max(A) / sigma_min(A) is predicted
01506 *        to be greater than the overflow threshold. This is because the
01507 *        a posteriori computation of the singular vectors assumes robust
01508 *        implementation of BLAS and some LAPACK procedures, capable of working
01509 *        in presence of extreme values. Since that is not always the case, ...
01510 *
01511          DO 7968 p = 1, NR
01512             CALL SCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
01513  7968    CONTINUE
01514 *
01515          IF ( L2PERT ) THEN
01516             XSC = SQRT(SMALL/EPSLN)
01517             DO 5969 q = 1, NR
01518                TEMP1 = XSC*ABS( V(q,q) )
01519                DO 5968 p = 1, N
01520                   IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
01521      &                .OR. ( p .LT. q ) )
01522      &                V(p,q) = SIGN( TEMP1, V(p,q) )
01523                   IF ( p. LT. q ) V(p,q) = - V(p,q)
01524  5968          CONTINUE
01525  5969       CONTINUE
01526          ELSE
01527             CALL SLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01528          END IF
01529 
01530          CALL SGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
01531      &        LWORK-2*N, IERR )
01532          CALL SLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N )
01533 *
01534          DO 7969 p = 1, NR
01535             CALL SCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
01536  7969    CONTINUE
01537 
01538          IF ( L2PERT ) THEN
01539             XSC = SQRT(SMALL/EPSLN)
01540             DO 9970 q = 2, NR
01541                DO 9971 p = 1, q - 1
01542                   TEMP1 = XSC * AMIN1(ABS(U(p,p)),ABS(U(q,q)))
01543                   U(p,q) = - SIGN( TEMP1, U(q,p) )
01544  9971          CONTINUE
01545  9970       CONTINUE
01546          ELSE
01547             CALL SLASET('U', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
01548          END IF
01549 
01550          CALL SGESVJ( 'L', 'U', 'V', NR, NR, U, LDU, SVA,
01551      &        N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO )
01552          SCALEM  = WORK(2*N+N*NR+1)
01553          NUMRANK = NINT(WORK(2*N+N*NR+2))
01554 
01555          IF ( NR .LT. N ) THEN
01556             CALL SLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
01557             CALL SLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
01558             CALL SLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
01559          END IF
01560 
01561          CALL SORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
01562      &        V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
01563 *
01564 *           Permute the rows of V using the (column) permutation from the
01565 *           first QRF. Also, scale the columns to make them unit in
01566 *           Euclidean norm. This applies to all cases.
01567 *
01568             TEMP1 = SQRT(FLOAT(N)) * EPSLN
01569             DO 7972 q = 1, N
01570                DO 8972 p = 1, N
01571                   WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
01572  8972          CONTINUE
01573                DO 8973 p = 1, N
01574                   V(p,q) = WORK(2*N+N*NR+NR+p)
01575  8973          CONTINUE
01576                XSC = ONE / SNRM2( N, V(1,q), 1 )
01577                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01578      &           CALL SSCAL( N, XSC, V(1,q), 1 )
01579  7972       CONTINUE
01580 *
01581 *           At this moment, V contains the right singular vectors of A.
01582 *           Next, assemble the left singular vector matrix U (M x N).
01583 *
01584          IF ( NR .LT. M ) THEN
01585             CALL SLASET( 'A',  M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
01586             IF ( NR .LT. N1 ) THEN
01587                CALL SLASET( 'A',NR,  N1-NR, ZERO, ZERO,  U(1,NR+1),LDU )
01588                CALL SLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
01589             END IF
01590          END IF
01591 *
01592          CALL SORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
01593      &        LDU, WORK(N+1), LWORK-N, IERR )
01594 *
01595             IF ( ROWPIV )
01596      &         CALL SLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
01597 *
01598 *
01599          END IF
01600          IF ( TRANSP ) THEN
01601 *           .. swap U and V because the procedure worked on A^t
01602             DO 6974 p = 1, N
01603                CALL SSWAP( N, U(1,p), 1, V(1,p), 1 )
01604  6974       CONTINUE
01605          END IF
01606 *
01607       END IF
01608 *     end of the full SVD
01609 *
01610 *     Undo scaling, if necessary (and possible)
01611 *
01612       IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
01613          CALL SLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
01614          USCAL1 = ONE
01615          USCAL2 = ONE
01616       END IF
01617 *
01618       IF ( NR .LT. N ) THEN
01619          DO 3004 p = NR+1, N
01620             SVA(p) = ZERO
01621  3004    CONTINUE
01622       END IF
01623 *
01624       WORK(1) = USCAL2 * SCALEM
01625       WORK(2) = USCAL1
01626       IF ( ERREST ) WORK(3) = SCONDA
01627       IF ( LSVEC .AND. RSVEC ) THEN
01628          WORK(4) = CONDR1
01629          WORK(5) = CONDR2
01630       END IF
01631       IF ( L2TRAN ) THEN
01632          WORK(6) = ENTRA
01633          WORK(7) = ENTRAT
01634       END IF
01635 *
01636       IWORK(1) = NR
01637       IWORK(2) = NUMRANK
01638       IWORK(3) = WARNING
01639 *
01640       RETURN
01641 *     ..
01642 *     .. END OF SGEJSV
01643 *     ..
01644       END
01645 *