--- rpl/lapack/lapack/dgesvj.f 2010/08/07 13:22:14 1.2 +++ rpl/lapack/lapack/dgesvj.f 2020/05/21 21:45:57 1.20 @@ -1,283 +1,382 @@ - SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, - + LDV, WORK, LWORK, INFO ) +*> \brief \b DGESVJ +* +* =========== DOCUMENTATION =========== * -* -- LAPACK routine (version 3.2.2) -- +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ * -* -- Contributed by Zlatko Drmac of the University of Zagreb and -- -* -- Kresimir Veselic of the Fernuniversitaet Hagen -- -* -- June 2010 -- +*> \htmlonly +*> Download DGESVJ + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, +* LDV, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LDA, LDV, LWORK, M, MV, N +* CHARACTER*1 JOBA, JOBU, JOBV +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ), +* $ WORK( LWORK ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGESVJ 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^t, [++] = [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. +*> DGESVJ can sometimes compute tiny singular values and their singular vectors much +*> more accurately than other SVD routines, see below under Further Details. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBA +*> \verbatim +*> JOBA is CHARACTER*1 +*> Specifies the structure of A. +*> = 'L': The input matrix A is lower triangular; +*> = 'U': The input matrix A is upper triangular; +*> = 'G': The input matrix A is general M-by-N matrix, M >= N. +*> \endverbatim +*> +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> Specifies whether to compute the left singular vectors +*> (columns of U): +*> = 'U': The left singular vectors corresponding to the nonzero +*> singular values are computed and returned in the leading +*> columns of A. See more details in the description of A. +*> The default numerical orthogonality threshold is set to +*> approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E'). +*> = 'C': Analogous to JOBU='U', except that user can control the +*> level of numerical orthogonality of the computed left +*> singular vectors. TOL can be set to TOL = CTOL*EPS, where +*> CTOL is given on input in the array WORK. +*> No CTOL smaller than ONE is allowed. CTOL greater +*> than 1 / EPS is meaningless. The option 'C' +*> can be used if M*EPS is satisfactory orthogonality +*> of the computed left singular vectors, so CTOL=M could +*> save few sweeps of Jacobi rotations. +*> See the descriptions of A and WORK(1). +*> = 'N': The matrix U is not computed. However, see the +*> description of A. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> 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 +*> 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 +*> referenced +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the input matrix A. +*> M >= N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit : +*> 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 +*> above the underflow threshold DLAMCH('S'). The singular +*> vectors corresponding to underflowed or zero singular +*> values are not computed. The value of RANKA is returned +*> in the array WORK as RANKA=NINT(WORK(2)). Also see the +*> descriptions of SVA and WORK. The computed columns of U +*> are mutually numerically orthogonal up to approximately +*> TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), +*> see the description of JOBU. +*> If INFO > 0 : +*> the procedure DGESVJ did not converge in the given number +*> of iterations (sweeps). In that case, the computed +*> columns of U may not be orthogonal up to TOL. The output +*> U (stored in A), SIGMA (given by the computed singular +*> 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 = '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 +*> orthogonality of U is not an issue and iterations are +*> stopped when the columns of the iterated matrix are +*> numerically orthogonal up to approximately M*EPS. Thus, +*> on exit, A contains the columns of U scaled with the +*> corresponding singular values. +*> If INFO > 0 : +*> the procedure DGESVJ did not converge in the given number +*> of iterations (sweeps). +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[out] SVA +*> \verbatim +*> SVA is DOUBLE PRECISION array, dimension (N) +*> On exit : +*> If INFO = 0 : +*> depending on the value SCALE = WORK(1), we have: +*> 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. +*> If SCALE .NE. ONE : +*> The singular values of A are SCALE*SVA(1:N), and this +*> factored representation is due to the fact that some of the +*> singular values of A might underflow or overflow. +*> If INFO > 0 : +*> the procedure DGESVJ did not converge in the given number of +*> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate. +*> \endverbatim +*> +*> \param[in] MV +*> \verbatim +*> MV is INTEGER +*> If JOBV = 'A', then the product of Jacobi rotations in DGESVJ +*> is applied to the first MV rows of V. See the description of JOBV. +*> \endverbatim +*> +*> \param[in,out] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension (LDV,N) +*> If JOBV = 'V', then V contains on exit the N-by-N matrix of +*> the right singular vectors; +*> If JOBV = 'A', then V contains the product of the computed right +*> singular vector matrix and the initial matrix in +*> the array V. +*> If JOBV = 'N', then V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V, LDV >= 1. +*> If JOBV = 'V', then LDV >= max(1,N). +*> If JOBV = 'A', then LDV >= max(1,MV) . +*> \endverbatim +*> +*> \param[in,out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (LWORK) +*> On entry : +*> If JOBU = 'C' : +*> WORK(1) = CTOL, where CTOL defines the threshold for convergence. +*> The process stops if all columns of A are mutually +*> orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). +*> It is required that CTOL >= ONE, i.e. it is not +*> allowed to force the routine to obtain orthogonality +*> below EPS. +*> On exit : +*> WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) +*> are the computed singular values of A. +*> (See description of SVA().) +*> WORK(2) = NINT(WORK(2)) is the number of the computed nonzero +*> singular values. +*> WORK(3) = NINT(WORK(3)) is the number of the computed singular +*> values that are larger than the underflow threshold. +*> WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi +*> rotations needed for numerical convergence. +*> WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. +*> This is useful information in cases when DGESVJ did +*> not converge, as it can be used to estimate whether +*> the output is 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. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> length of WORK, WORK >= MAX(6,M+N) +*> \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: DGESVJ did not converge in the maximal allowed number (30) +*> of sweeps. The output may still be useful. See the +*> description of WORK. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date June 2017 +* +*> \ingroup doubleGEcomputational +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane +*> rotations. The rotations are implemented as fast scaled rotations of +*> Anda and Park [1]. In the case of underflow of the Jacobi angle, a +*> modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses +*> column interchanges of de Rijk [2]. The relative accuracy of the computed +*> singular values and the accuracy of the computed singular vectors (in +*> angle metric) is as guaranteed by the theory of Demmel and Veselic [3]. +*> The condition number that determines the accuracy in the full rank case +*> is essentially min_{D=diag} kappa(A*D), where kappa(.) is the +*> spectral condition number. The best performance of this Jacobi SVD +*> procedure is achieved if used in an accelerated version of Drmac and +*> Veselic [5,6], and it is the kernel routine in the SIGMA library [7]. +*> Some tunning parameters (marked with [TP]) are available for the +*> implementer. +*> The computational range for the nonzero singular values is the machine +*> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even +*> denormalized singular values can be computed with the corresponding +*> gradual loss of accurate digits. +*> \endverbatim +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> +*> ============ +*> +*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) +*> \endverbatim +* +*> \par References: +* ================ +*> +*> \verbatim +*> +*> [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling. +*> SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174. +*> [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the +*> singular value decomposition on a vector computer. +*> SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. +*> [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. +*> [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular +*> value computation in floating point arithmetic. +*> SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. +*> [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. +*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. +*> LAPACK Working note 169. +*> [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. +*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. +*> LAPACK Working note 170. +*> [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, +*> QSVD, (H,K)-SVD computations. +*> Department of Mathematics, University of Zagreb, 2008. +*> \endverbatim +* +*> \par Bugs, examples and comments: +* ================================= +*> +*> \verbatim +*> =========================== +*> Please report all bugs and send interesting test examples and comments to +*> drmac@math.hr. Thank you. +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, + $ LDV, WORK, LWORK, INFO ) * +* -- LAPACK computational routine (version 3.7.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* June 2017 * -* This routine is also part of SIGMA (version 1.23, October 23. 2008.) -* SIGMA is a library of algorithms for highly accurate algorithms for -* computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the -* eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0. -* - IMPLICIT NONE -* .. * .. Scalar Arguments .. INTEGER INFO, LDA, LDV, LWORK, M, MV, N CHARACTER*1 JOBA, JOBU, JOBV * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ), - + WORK( LWORK ) + $ WORK( LWORK ) * .. * -* Purpose -* ======= -* -* DGESVJ 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^t, [++] = [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. -* -* Further Details -* ~~~~~~~~~~~~~~~ -* The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane -* rotations. The rotations are implemented as fast scaled rotations of -* Anda and Park [1]. In the case of underflow of the Jacobi angle, a -* modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses -* column interchanges of de Rijk [2]. The relative accuracy of the computed -* singular values and the accuracy of the computed singular vectors (in -* angle metric) is as guaranteed by the theory of Demmel and Veselic [3]. -* The condition number that determines the accuracy in the full rank case -* is essentially min_{D=diag} kappa(A*D), where kappa(.) is the -* spectral condition number. The best performance of this Jacobi SVD -* procedure is achieved if used in an accelerated version of Drmac and -* Veselic [5,6], and it is the kernel routine in the SIGMA library [7]. -* Some tunning parameters (marked with [TP]) are available for the -* implementer. -* The computational range for the nonzero singular values is the machine -* number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even -* denormalized singular values can be computed with the corresponding -* gradual loss of accurate digits. -* -* Contributors -* ~~~~~~~~~~~~ -* Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) -* -* References -* ~~~~~~~~~~ -* [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling. -* SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174. -* [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the -* singular value decomposition on a vector computer. -* SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. -* [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. -* [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular -* value computation in floating point arithmetic. -* SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. -* [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. -* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. -* LAPACK Working note 169. -* [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. -* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. -* LAPACK Working note 170. -* [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, -* QSVD, (H,K)-SVD computations. -* Department of Mathematics, University of Zagreb, 2008. -* -* Bugs, Examples and Comments -* ~~~~~~~~~~~~~~~~~~~~~~~~~~~ -* Please report all bugs and send interesting test examples and comments to -* drmac@math.hr. Thank you. -* -* Arguments -* ========= -* -* JOBA (input) CHARACTER* 1 -* Specifies the structure of A. -* = 'L': The input matrix A is lower triangular; -* = 'U': The input matrix A is upper triangular; -* = 'G': The input matrix A is general M-by-N matrix, M >= N. -* -* JOBU (input) CHARACTER*1 -* Specifies whether to compute the left singular vectors -* (columns of U): -* = 'U': The left singular vectors corresponding to the nonzero -* singular values are computed and returned in the leading -* columns of A. See more details in the description of A. -* The default numerical orthogonality threshold is set to -* approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E'). -* = 'C': Analogous to JOBU='U', except that user can control the -* level of numerical orthogonality of the computed left -* singular vectors. TOL can be set to TOL = CTOL*EPS, where -* CTOL is given on input in the array WORK. -* No CTOL smaller than ONE is allowed. CTOL greater -* than 1 / EPS is meaningless. The option 'C' -* can be used if M*EPS is satisfactory orthogonality -* of the computed left singular vectors, so CTOL=M could -* save few sweeps of Jacobi rotations. -* See the descriptions of A and WORK(1). -* = 'N': The matrix U is not computed. However, see the -* description of A. -* -* JOBV (input) 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 -* 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 -* referenced -* -* M (input) INTEGER -* The number of rows of the input matrix A. M >= 0. -* -* N (input) INTEGER -* The number of columns of the input matrix A. -* M >= N >= 0. -* -* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) -* On entry, the M-by-N matrix A. -* On exit : -* If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' : -* If INFO .EQ. 0 : -* 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 -* above the underflow threshold DLAMCH('S'). The singular -* vectors corresponding to underflowed or zero singular -* values are not computed. The value of RANKA is returned -* in the array WORK as RANKA=NINT(WORK(2)). Also see the -* descriptions of SVA and WORK. The computed columns of U -* are mutually numerically orthogonal up to approximately -* TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'), -* see the description of JOBU. -* If INFO .GT. 0 : -* the procedure DGESVJ did not converge in the given number -* of iterations (sweeps). In that case, the computed -* columns of U may not be orthogonal up to TOL. The output -* U (stored in A), SIGMA (given by the computed singular -* 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 : -* 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 -* orthogonality of U is not an issue and iterations are -* stopped when the columns of the iterated matrix are -* 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 : -* the procedure DGESVJ did not converge in the given number -* of iterations (sweeps). -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,M). -* -* SVA (workspace/output) DOUBLE PRECISION array, dimension (N) -* On exit : -* If INFO .EQ. 0 : -* depending on the value SCALE = WORK(1), we have: -* If SCALE .EQ. 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. -* If SCALE .NE. ONE : -* The singular values of A are SCALE*SVA(1:N), and this -* factored representation is due to the fact that some of the -* singular values of A might underflow or overflow. -* If INFO .GT. 0 : -* the procedure DGESVJ did not converge in the given number of -* iterations (sweeps) and SCALE*SVA(1:N) may not be accurate. -* -* MV (input) INTEGER -* If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ -* is applied to the first MV rows of V. See the description of JOBV. -* -* V (input/output) DOUBLE PRECISION array, dimension (LDV,N) -* If JOBV = 'V', then V contains on exit the N-by-N matrix of -* the right singular vectors; -* If JOBV = 'A', then V contains the product of the computed right -* singular vector matrix and the initial matrix in -* the array V. -* If JOBV = 'N', then V is not referenced. -* -* LDV (input) 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) . -* -* WORK (input/workspace/output) DOUBLE PRECISION array, dimension max(4,M+N). -* On entry : -* If JOBU .EQ. 'C' : -* WORK(1) = CTOL, where CTOL defines the threshold for convergence. -* The process stops if all columns of A are mutually -* orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). -* It is required that CTOL >= ONE, i.e. it is not -* allowed to force the routine to obtain orthogonality -* below EPSILON. -* On exit : -* WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) -* are the computed singular values of A. -* (See description of SVA().) -* WORK(2) = NINT(WORK(2)) is the number of the computed nonzero -* singular values. -* WORK(3) = NINT(WORK(3)) is the number of the computed singular -* values that are larger than the underflow threshold. -* WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi -* rotations needed for numerical convergence. -* WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. -* This is useful information in cases when DGESVJ did -* not converge, as it can be used to estimate whether -* the output is stil useful and for post festum analysis. -* 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. -* -* LWORK (input) INTEGER -* length of WORK, WORK >= MAX(6,M+N) -* -* INFO (output) INTEGER -* = 0 : successful exit. -* < 0 : if INFO = -i, then the i-th argument had an illegal value -* > 0 : DGESVJ did not converge in the maximal allowed number (30) -* of sweeps. The output may still be useful. See the -* description of WORK. -* * ===================================================================== * * .. Local Parameters .. - DOUBLE PRECISION ZERO, HALF, ONE, TWO - PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, - + TWO = 2.0D0 ) + DOUBLE PRECISION ZERO, HALF, ONE + PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) INTEGER NSWEEP PARAMETER ( NSWEEP = 30 ) * .. * .. Local Scalars .. DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG, - + BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ, - + MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, - + SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA, - + THSIGN, TOL + $ BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ, + $ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, + $ SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, + $ THSIGN, TOL INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, - + ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, - + N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, - + SWBAND + $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, + $ N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, + $ SWBAND LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK, - + RSVEC, UCTOL, UPPER + $ RSVEC, UCTOL, UPPER * .. * .. Local Arrays .. DOUBLE PRECISION FASTR( 5 ) * .. * .. Intrinsic Functions .. - INTRINSIC DABS, DMAX1, DMIN1, DBLE, MIN0, DSIGN, DSQRT + INTRINSIC DABS, MAX, MIN, DBLE, DSIGN, DSQRT * .. * .. External Functions .. * .. @@ -327,11 +426,11 @@ ELSE IF( MV.LT.0 ) THEN INFO = -9 ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR. - + ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN + $ ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN INFO = -11 ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN INFO = -12 - ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN + ELSE IF( LWORK.LT.MAX( M+N, 6 ) ) THEN INFO = -13 ELSE INFO = 0 @@ -368,22 +467,22 @@ * ... and the machine dependent parameters are *[!] (Make sure that DLAMCH() works properly on the target machine.) * - EPSILON = DLAMCH( 'Epsilon' ) - ROOTEPS = DSQRT( EPSILON ) + EPSLN = DLAMCH( 'Epsilon' ) + ROOTEPS = DSQRT( EPSLN ) SFMIN = DLAMCH( 'SafeMinimum' ) ROOTSFMIN = DSQRT( SFMIN ) - SMALL = SFMIN / EPSILON + SMALL = SFMIN / EPSLN BIG = DLAMCH( 'Overflow' ) * BIG = ONE / SFMIN ROOTBIG = ONE / ROOTSFMIN LARGE = BIG / DSQRT( DBLE( M*N ) ) BIGTHETA = ONE / ROOTEPS * - TOL = CTOL*EPSILON + TOL = CTOL*EPSLN ROOTTOL = DSQRT( TOL ) * - IF( DBLE( M )*EPSILON.GE.ONE ) THEN - INFO = -5 + IF( DBLE( M )*EPSLN.GE.ONE ) THEN + INFO = -4 CALL XERBLA( 'DGESVJ', -INFO ) RETURN END IF @@ -407,7 +506,7 @@ * DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries * in A are detected, the procedure returns with INFO=-6. * - SCALE = ONE / DSQRT( DBLE( M )*DBLE( N ) ) + SKL= ONE / DSQRT( DBLE( M )*DBLE( N ) ) NOSCALE = .TRUE. GOSCALE = .TRUE. * @@ -415,7 +514,7 @@ * the input matrix is M-by-N lower triangular (trapezoidal) DO 1874 p = 1, N AAPP = ZERO - AAQQ = ZERO + AAQQ = ONE CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ ) IF( AAPP.GT.BIG ) THEN INFO = -6 @@ -427,11 +526,11 @@ SVA( p ) = AAPP*AAQQ ELSE NOSCALE = .FALSE. - SVA( p ) = AAPP*( AAQQ*SCALE ) + SVA( p ) = AAPP*( AAQQ*SKL) IF( GOSCALE ) THEN GOSCALE = .FALSE. DO 1873 q = 1, p - 1 - SVA( q ) = SVA( q )*SCALE + SVA( q ) = SVA( q )*SKL 1873 CONTINUE END IF END IF @@ -440,7 +539,7 @@ * the input matrix is M-by-N upper triangular (trapezoidal) DO 2874 p = 1, N AAPP = ZERO - AAQQ = ZERO + AAQQ = ONE CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ ) IF( AAPP.GT.BIG ) THEN INFO = -6 @@ -452,11 +551,11 @@ SVA( p ) = AAPP*AAQQ ELSE NOSCALE = .FALSE. - SVA( p ) = AAPP*( AAQQ*SCALE ) + SVA( p ) = AAPP*( AAQQ*SKL) IF( GOSCALE ) THEN GOSCALE = .FALSE. DO 2873 q = 1, p - 1 - SVA( q ) = SVA( q )*SCALE + SVA( q ) = SVA( q )*SKL 2873 CONTINUE END IF END IF @@ -465,7 +564,7 @@ * the input matrix is M-by-N general dense DO 3874 p = 1, N AAPP = ZERO - AAQQ = ZERO + AAQQ = ONE CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ ) IF( AAPP.GT.BIG ) THEN INFO = -6 @@ -477,18 +576,18 @@ SVA( p ) = AAPP*AAQQ ELSE NOSCALE = .FALSE. - SVA( p ) = AAPP*( AAQQ*SCALE ) + SVA( p ) = AAPP*( AAQQ*SKL) IF( GOSCALE ) THEN GOSCALE = .FALSE. DO 3873 q = 1, p - 1 - SVA( q ) = SVA( q )*SCALE + SVA( q ) = SVA( q )*SKL 3873 CONTINUE END IF END IF 3874 CONTINUE END IF * - IF( NOSCALE )SCALE = ONE + IF( NOSCALE )SKL= ONE * * Move the smaller part of the spectrum from the underflow threshold *(!) Start by determining the position of the nonzero entries of the @@ -497,8 +596,8 @@ AAPP = ZERO AAQQ = BIG DO 4781 p = 1, N - IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) ) - AAPP = DMAX1( AAPP, SVA( p ) ) + IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) ) + AAPP = MAX( AAPP, SVA( p ) ) 4781 CONTINUE * * #:) Quick return for zero matrix @@ -517,9 +616,9 @@ * #:) Quick return for one-column matrix * IF( N.EQ.1 ) THEN - IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1, - + A( 1, 1 ), LDA, IERR ) - WORK( 1 ) = ONE / SCALE + IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1, + $ A( 1, 1 ), LDA, IERR ) + WORK( 1 ) = ONE / SKL IF( SVA( 1 ).GE.SFMIN ) THEN WORK( 2 ) = ONE ELSE @@ -535,23 +634,23 @@ * Protect small singular values from underflow, and try to * avoid underflows/overflows in computing Jacobi rotations. * - SN = DSQRT( SFMIN / EPSILON ) + SN = DSQRT( SFMIN / EPSLN ) TEMP1 = DSQRT( BIG / DBLE( N ) ) IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. - + ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN - TEMP1 = DMIN1( BIG, TEMP1 / AAPP ) + $ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN + TEMP1 = MIN( BIG, TEMP1 / AAPP ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN - TEMP1 = DMIN1( SN / AAQQ, BIG / ( AAPP*DSQRT( DBLE( N ) ) ) ) + TEMP1 = MIN( SN / AAQQ, BIG / ( AAPP*DSQRT( DBLE( N ) ) ) ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN - TEMP1 = DMAX1( SN / AAQQ, TEMP1 / AAPP ) + TEMP1 = MAX( SN / AAQQ, TEMP1 / AAPP ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN - TEMP1 = DMIN1( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) ) + TEMP1 = MIN( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) ) * AAQQ = AAQQ*TEMP1 * AAPP = AAPP*TEMP1 ELSE @@ -563,10 +662,10 @@ IF( TEMP1.NE.ONE ) THEN CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR ) END IF - SCALE = TEMP1*SCALE - IF( SCALE.NE.ONE ) THEN - CALL DLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR ) - SCALE = ONE / SCALE + SKL= TEMP1*SKL + IF( SKL.NE.ONE ) THEN + CALL DLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR ) + SKL= ONE / SKL END IF * * Row-cyclic Jacobi SVD algorithm with column pivoting @@ -592,7 +691,7 @@ * The boundaries are determined dynamically, based on the number of * pivots above a threshold. * - KBL = MIN0( 8, N ) + KBL = MIN( 8, N ) *[TP] KBL is a tuning parameter that defines the tile size in the * tiling of the p-q loops of pivot pairs. In general, an optimal * value of KBL depends on the matrix dimensions and on the @@ -604,7 +703,7 @@ BLSKIP = KBL**2 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL. * - ROWSKIP = MIN0( 5, KBL ) + ROWSKIP = MIN( 5, KBL ) *[TP] ROWSKIP is a tuning parameter. * LKAHEAD = 1 @@ -615,7 +714,7 @@ * invokes cubic convergence. Big part of this cycle is done inside * canonical subspaces of dimensions less than M. * - IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN + IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX( 64, 4*KBL ) ) ) THEN *[TP] The number of partition levels and the actual partition are * tuning parameters. N4 = N / 4 @@ -638,54 +737,54 @@ * [+ + x x] [x x]. [x x] * CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA, - + WORK( N34+1 ), SVA( N34+1 ), MVL, - + V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL, - + 2, WORK( N+1 ), LWORK-N, IERR ) + $ WORK( N34+1 ), SVA( N34+1 ), MVL, + $ V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL, + $ 2, WORK( N+1 ), LWORK-N, IERR ) * CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA, - + WORK( N2+1 ), SVA( N2+1 ), MVL, - + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2, - + WORK( N+1 ), LWORK-N, IERR ) + $ WORK( N2+1 ), SVA( N2+1 ), MVL, + $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2, + $ WORK( N+1 ), LWORK-N, IERR ) * CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA, - + WORK( N2+1 ), SVA( N2+1 ), MVL, - + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1, - + WORK( N+1 ), LWORK-N, IERR ) + $ WORK( N2+1 ), SVA( N2+1 ), MVL, + $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, + $ WORK( N+1 ), LWORK-N, IERR ) * CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA, - + WORK( N4+1 ), SVA( N4+1 ), MVL, - + V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1, - + WORK( N+1 ), LWORK-N, IERR ) + $ WORK( N4+1 ), SVA( N4+1 ), MVL, + $ V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1, + $ WORK( N+1 ), LWORK-N, IERR ) * CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV, - + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N, - + IERR ) + $ EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N, + $ IERR ) * CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V, - + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ), - + LWORK-N, IERR ) + $ LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ), + $ LWORK-N, IERR ) * * ELSE IF( UPPER ) THEN * * CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV, - + EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N, - + IERR ) + $ EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N, + $ IERR ) * CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ), - + SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV, - + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N, - + IERR ) + $ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV, + $ EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N, + $ IERR ) * CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V, - + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ), - + LWORK-N, IERR ) + $ LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ), + $ LWORK-N, IERR ) * CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA, - + WORK( N2+1 ), SVA( N2+1 ), MVL, - + V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1, - + WORK( N+1 ), LWORK-N, IERR ) + $ WORK( N2+1 ), SVA( N2+1 ), MVL, + $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, + $ WORK( N+1 ), LWORK-N, IERR ) END IF * @@ -713,11 +812,11 @@ * igl = ( ibr-1 )*KBL + 1 * - DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr ) + DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr ) * igl = igl + ir1*KBL * - DO 2001 p = igl, MIN0( igl+KBL-1, N-1 ) + DO 2001 p = igl, MIN( igl+KBL-1, N-1 ) * * .. de Rijk's pivoting * @@ -725,7 +824,7 @@ IF( p.NE.q ) THEN CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, - + V( 1, q ), 1 ) + $ V( 1, q ), 1 ) TEMP1 = SVA( p ) SVA( p ) = SVA( q ) SVA( q ) = TEMP1 @@ -749,11 +848,11 @@ * below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)". * IF( ( SVA( p ).LT.ROOTBIG ) .AND. - + ( SVA( p ).GT.ROOTSFMIN ) ) THEN + $ ( SVA( p ).GT.ROOTSFMIN ) ) THEN SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p ) ELSE TEMP1 = ZERO - AAPP = ZERO + AAPP = ONE CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p ) END IF @@ -766,7 +865,7 @@ * PSKIPPED = 0 * - DO 2002 q = p + 1, MIN0( igl+KBL-1, N ) + DO 2002 q = p + 1, MIN( igl+KBL-1, N ) * AAQQ = SVA( q ) * @@ -777,35 +876,35 @@ ROTOK = ( SMALL*AAPP ).LE.AAQQ IF( AAPP.LT.( BIG / AAQQ ) ) THEN AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, - + q ), 1 )*WORK( p )*WORK( q ) / - + AAQQ ) / AAPP + $ q ), 1 )*WORK( p )*WORK( q ) / + $ AAQQ ) / AAPP ELSE CALL DCOPY( M, A( 1, p ), 1, - + WORK( N+1 ), 1 ) + $ WORK( N+1 ), 1 ) CALL DLASCL( 'G', 0, 0, AAPP, - + WORK( p ), M, 1, - + WORK( N+1 ), LDA, IERR ) + $ WORK( p ), M, 1, + $ WORK( N+1 ), LDA, IERR ) AAPQ = DDOT( M, WORK( N+1 ), 1, - + A( 1, q ), 1 )*WORK( q ) / AAQQ + $ A( 1, q ), 1 )*WORK( q ) / AAQQ END IF ELSE ROTOK = AAPP.LE.( AAQQ / SMALL ) IF( AAPP.GT.( SMALL / AAQQ ) ) THEN AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, - + q ), 1 )*WORK( p )*WORK( q ) / - + AAQQ ) / AAPP + $ q ), 1 )*WORK( p )*WORK( q ) / + $ AAQQ ) / AAPP ELSE CALL DCOPY( M, A( 1, q ), 1, - + WORK( N+1 ), 1 ) + $ WORK( N+1 ), 1 ) CALL DLASCL( 'G', 0, 0, AAQQ, - + WORK( q ), M, 1, - + WORK( N+1 ), LDA, IERR ) + $ WORK( q ), M, 1, + $ WORK( N+1 ), LDA, IERR ) AAPQ = DDOT( M, WORK( N+1 ), 1, - + A( 1, p ), 1 )*WORK( p ) / AAPP + $ A( 1, p ), 1 )*WORK( p ) / AAPP END IF END IF * - MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) ) + MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) ) * * TO rotate or NOT to rotate, THAT is the question ... * @@ -824,26 +923,25 @@ * AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ - THETA = -HALF*DABS( AQOAP-APOAQ ) / - + AAPQ + THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ * IF( DABS( THETA ).GT.BIGTHETA ) THEN * T = HALF / THETA FASTR( 3 ) = T*WORK( p ) / WORK( q ) FASTR( 4 ) = -T*WORK( q ) / - + WORK( p ) + $ WORK( p ) CALL DROTM( M, A( 1, p ), 1, - + A( 1, q ), 1, FASTR ) + $ A( 1, q ), 1, FASTR ) IF( RSVEC )CALL DROTM( MVL, - + V( 1, p ), 1, - + V( 1, q ), 1, - + FASTR ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, - + ONE+T*APOAQ*AAPQ ) ) - AAPP = AAPP*DSQRT( ONE-T*AQOAP* - + AAPQ ) - MXSINJ = DMAX1( MXSINJ, DABS( T ) ) + $ V( 1, p ), 1, + $ V( 1, q ), 1, + $ FASTR ) + SVA( q ) = AAQQ*DSQRT( MAX( ZERO, + $ ONE+T*APOAQ*AAPQ ) ) + AAPP = AAPP*DSQRT( MAX( ZERO, + $ ONE-T*AQOAP*AAPQ ) ) + MXSINJ = MAX( MXSINJ, DABS( T ) ) * ELSE * @@ -851,15 +949,15 @@ * THSIGN = -DSIGN( ONE, AAPQ ) T = ONE / ( THETA+THSIGN* - + DSQRT( ONE+THETA*THETA ) ) + $ DSQRT( ONE+THETA*THETA ) ) CS = DSQRT( ONE / ( ONE+T*T ) ) SN = T*CS * - MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, - + ONE+T*APOAQ*AAPQ ) ) - AAPP = AAPP*DSQRT( DMAX1( ZERO, - + ONE-T*AQOAP*AAPQ ) ) + MXSINJ = MAX( MXSINJ, DABS( SN ) ) + SVA( q ) = AAQQ*DSQRT( MAX( ZERO, + $ ONE+T*APOAQ*AAPQ ) ) + AAPP = AAPP*DSQRT( MAX( ZERO, + $ ONE-T*AQOAP*AAPQ ) ) * APOAQ = WORK( p ) / WORK( q ) AQOAP = WORK( q ) / WORK( p ) @@ -870,88 +968,88 @@ WORK( p ) = WORK( p )*CS WORK( q ) = WORK( q )*CS CALL DROTM( M, A( 1, p ), 1, - + A( 1, q ), 1, - + FASTR ) + $ A( 1, q ), 1, + $ FASTR ) IF( RSVEC )CALL DROTM( MVL, - + V( 1, p ), 1, V( 1, q ), - + 1, FASTR ) + $ V( 1, p ), 1, V( 1, q ), + $ 1, FASTR ) ELSE CALL DAXPY( M, -T*AQOAP, - + A( 1, q ), 1, - + A( 1, p ), 1 ) + $ A( 1, q ), 1, + $ A( 1, p ), 1 ) CALL DAXPY( M, CS*SN*APOAQ, - + A( 1, p ), 1, - + A( 1, q ), 1 ) + $ A( 1, p ), 1, + $ A( 1, q ), 1 ) WORK( p ) = WORK( p )*CS WORK( q ) = WORK( q ) / CS IF( RSVEC ) THEN CALL DAXPY( MVL, -T*AQOAP, - + V( 1, q ), 1, - + V( 1, p ), 1 ) + $ V( 1, q ), 1, + $ V( 1, p ), 1 ) CALL DAXPY( MVL, - + CS*SN*APOAQ, - + V( 1, p ), 1, - + V( 1, q ), 1 ) + $ CS*SN*APOAQ, + $ V( 1, p ), 1, + $ V( 1, q ), 1 ) END IF END IF ELSE IF( WORK( q ).GE.ONE ) THEN CALL DAXPY( M, T*APOAQ, - + A( 1, p ), 1, - + A( 1, q ), 1 ) + $ A( 1, p ), 1, + $ A( 1, q ), 1 ) CALL DAXPY( M, -CS*SN*AQOAP, - + A( 1, q ), 1, - + A( 1, p ), 1 ) + $ A( 1, q ), 1, + $ A( 1, p ), 1 ) WORK( p ) = WORK( p ) / CS WORK( q ) = WORK( q )*CS IF( RSVEC ) THEN CALL DAXPY( MVL, T*APOAQ, - + V( 1, p ), 1, - + V( 1, q ), 1 ) + $ V( 1, p ), 1, + $ V( 1, q ), 1 ) CALL DAXPY( MVL, - + -CS*SN*AQOAP, - + V( 1, q ), 1, - + V( 1, p ), 1 ) + $ -CS*SN*AQOAP, + $ V( 1, q ), 1, + $ V( 1, p ), 1 ) END IF ELSE IF( WORK( p ).GE.WORK( q ) ) - + THEN + $ THEN CALL DAXPY( M, -T*AQOAP, - + A( 1, q ), 1, - + A( 1, p ), 1 ) + $ A( 1, q ), 1, + $ A( 1, p ), 1 ) CALL DAXPY( M, CS*SN*APOAQ, - + A( 1, p ), 1, - + A( 1, q ), 1 ) + $ A( 1, p ), 1, + $ A( 1, q ), 1 ) WORK( p ) = WORK( p )*CS WORK( q ) = WORK( q ) / CS IF( RSVEC ) THEN CALL DAXPY( MVL, - + -T*AQOAP, - + V( 1, q ), 1, - + V( 1, p ), 1 ) + $ -T*AQOAP, + $ V( 1, q ), 1, + $ V( 1, p ), 1 ) CALL DAXPY( MVL, - + CS*SN*APOAQ, - + V( 1, p ), 1, - + V( 1, q ), 1 ) + $ CS*SN*APOAQ, + $ V( 1, p ), 1, + $ V( 1, q ), 1 ) END IF ELSE CALL DAXPY( M, T*APOAQ, - + A( 1, p ), 1, - + A( 1, q ), 1 ) + $ A( 1, p ), 1, + $ A( 1, q ), 1 ) CALL DAXPY( M, - + -CS*SN*AQOAP, - + A( 1, q ), 1, - + A( 1, p ), 1 ) + $ -CS*SN*AQOAP, + $ A( 1, q ), 1, + $ A( 1, p ), 1 ) WORK( p ) = WORK( p ) / CS WORK( q ) = WORK( q )*CS IF( RSVEC ) THEN CALL DAXPY( MVL, - + T*APOAQ, V( 1, p ), - + 1, V( 1, q ), 1 ) + $ T*APOAQ, V( 1, p ), + $ 1, V( 1, q ), 1 ) CALL DAXPY( MVL, - + -CS*SN*AQOAP, - + V( 1, q ), 1, - + V( 1, p ), 1 ) + $ -CS*SN*AQOAP, + $ V( 1, q ), 1, + $ V( 1, p ), 1 ) END IF END IF END IF @@ -961,20 +1059,20 @@ ELSE * .. have to use modified Gram-Schmidt like transformation CALL DCOPY( M, A( 1, p ), 1, - + WORK( N+1 ), 1 ) + $ WORK( N+1 ), 1 ) CALL DLASCL( 'G', 0, 0, AAPP, ONE, M, - + 1, WORK( N+1 ), LDA, - + IERR ) + $ 1, WORK( N+1 ), LDA, + $ IERR ) CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M, - + 1, A( 1, q ), LDA, IERR ) + $ 1, A( 1, q ), LDA, IERR ) TEMP1 = -AAPQ*WORK( p ) / WORK( q ) CALL DAXPY( M, TEMP1, WORK( N+1 ), 1, - + A( 1, q ), 1 ) + $ A( 1, q ), 1 ) CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M, - + 1, A( 1, q ), LDA, IERR ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, - + ONE-AAPQ*AAPQ ) ) - MXSINJ = DMAX1( MXSINJ, SFMIN ) + $ 1, A( 1, q ), LDA, IERR ) + SVA( q ) = AAQQ*DSQRT( MAX( ZERO, + $ ONE-AAPQ*AAPQ ) ) + MXSINJ = MAX( MXSINJ, SFMIN ) END IF * END IF ROTOK THEN ... ELSE * @@ -982,29 +1080,29 @@ * recompute SVA(q), SVA(p). * IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) - + THEN + $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. - + ( AAQQ.GT.ROOTSFMIN ) ) THEN + $ ( AAQQ.GT.ROOTSFMIN ) ) THEN SVA( q ) = DNRM2( M, A( 1, q ), 1 )* - + WORK( q ) + $ WORK( q ) ELSE T = ZERO - AAQQ = ZERO + AAQQ = ONE CALL DLASSQ( M, A( 1, q ), 1, T, - + AAQQ ) + $ AAQQ ) SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) END IF END IF IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN IF( ( AAPP.LT.ROOTBIG ) .AND. - + ( AAPP.GT.ROOTSFMIN ) ) THEN + $ ( AAPP.GT.ROOTSFMIN ) ) THEN AAPP = DNRM2( M, A( 1, p ), 1 )* - + WORK( p ) + $ WORK( p ) ELSE T = ZERO - AAPP = ZERO + AAPP = ONE CALL DLASSQ( M, A( 1, p ), 1, T, - + AAPP ) + $ AAPP ) AAPP = T*DSQRT( AAPP )*WORK( p ) END IF SVA( p ) = AAPP @@ -1023,7 +1121,7 @@ END IF * IF( ( i.LE.SWBAND ) .AND. - + ( PSKIPPED.GT.ROWSKIP ) ) THEN + $ ( PSKIPPED.GT.ROWSKIP ) ) THEN IF( ir1.EQ.0 )AAPP = -AAPP NOTROT = 0 GO TO 2103 @@ -1040,7 +1138,7 @@ ELSE SVA( p ) = AAPP IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) ) - + NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p + $ NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p END IF * 2001 CONTINUE @@ -1060,14 +1158,14 @@ * doing the block at ( ibr, jbc ) * IJBLSK = 0 - DO 2100 p = igl, MIN0( igl+KBL-1, N ) + DO 2100 p = igl, MIN( igl+KBL-1, N ) * AAPP = SVA( p ) IF( AAPP.GT.ZERO ) THEN * PSKIPPED = 0 * - DO 2200 q = jgl, MIN0( jgl+KBL-1, N ) + DO 2200 q = jgl, MIN( jgl+KBL-1, N ) * AAQQ = SVA( q ) IF( AAQQ.GT.ZERO ) THEN @@ -1085,16 +1183,16 @@ END IF IF( AAPP.LT.( BIG / AAQQ ) ) THEN AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, - + q ), 1 )*WORK( p )*WORK( q ) / - + AAQQ ) / AAPP + $ q ), 1 )*WORK( p )*WORK( q ) / + $ AAQQ ) / AAPP ELSE CALL DCOPY( M, A( 1, p ), 1, - + WORK( N+1 ), 1 ) + $ WORK( N+1 ), 1 ) CALL DLASCL( 'G', 0, 0, AAPP, - + WORK( p ), M, 1, - + WORK( N+1 ), LDA, IERR ) + $ WORK( p ), M, 1, + $ WORK( N+1 ), LDA, IERR ) AAPQ = DDOT( M, WORK( N+1 ), 1, - + A( 1, q ), 1 )*WORK( q ) / AAQQ + $ A( 1, q ), 1 )*WORK( q ) / AAQQ END IF ELSE IF( AAPP.GE.AAQQ ) THEN @@ -1104,20 +1202,20 @@ END IF IF( AAPP.GT.( SMALL / AAQQ ) ) THEN AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, - + q ), 1 )*WORK( p )*WORK( q ) / - + AAQQ ) / AAPP + $ q ), 1 )*WORK( p )*WORK( q ) / + $ AAQQ ) / AAPP ELSE CALL DCOPY( M, A( 1, q ), 1, - + WORK( N+1 ), 1 ) + $ WORK( N+1 ), 1 ) CALL DLASCL( 'G', 0, 0, AAQQ, - + WORK( q ), M, 1, - + WORK( N+1 ), LDA, IERR ) + $ WORK( q ), M, 1, + $ WORK( N+1 ), LDA, IERR ) AAPQ = DDOT( M, WORK( N+1 ), 1, - + A( 1, p ), 1 )*WORK( p ) / AAPP + $ A( 1, p ), 1 )*WORK( p ) / AAPP END IF END IF * - MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) ) + MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) ) * * TO rotate or NOT to rotate, THAT is the question ... * @@ -1131,26 +1229,25 @@ * AQOAP = AAQQ / AAPP APOAQ = AAPP / AAQQ - THETA = -HALF*DABS( AQOAP-APOAQ ) / - + AAPQ + THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ IF( AAQQ.GT.AAPP0 )THETA = -THETA * IF( DABS( THETA ).GT.BIGTHETA ) THEN T = HALF / THETA FASTR( 3 ) = T*WORK( p ) / WORK( q ) FASTR( 4 ) = -T*WORK( q ) / - + WORK( p ) + $ WORK( p ) CALL DROTM( M, A( 1, p ), 1, - + A( 1, q ), 1, FASTR ) + $ A( 1, q ), 1, FASTR ) IF( RSVEC )CALL DROTM( MVL, - + V( 1, p ), 1, - + V( 1, q ), 1, - + FASTR ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, - + ONE+T*APOAQ*AAPQ ) ) - AAPP = AAPP*DSQRT( DMAX1( ZERO, - + ONE-T*AQOAP*AAPQ ) ) - MXSINJ = DMAX1( MXSINJ, DABS( T ) ) + $ V( 1, p ), 1, + $ V( 1, q ), 1, + $ FASTR ) + SVA( q ) = AAQQ*DSQRT( MAX( ZERO, + $ ONE+T*APOAQ*AAPQ ) ) + AAPP = AAPP*DSQRT( MAX( ZERO, + $ ONE-T*AQOAP*AAPQ ) ) + MXSINJ = MAX( MXSINJ, DABS( T ) ) ELSE * * .. choose correct signum for THETA and rotate @@ -1158,14 +1255,14 @@ THSIGN = -DSIGN( ONE, AAPQ ) IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN T = ONE / ( THETA+THSIGN* - + DSQRT( ONE+THETA*THETA ) ) + $ DSQRT( ONE+THETA*THETA ) ) CS = DSQRT( ONE / ( ONE+T*T ) ) SN = T*CS - MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, - + ONE+T*APOAQ*AAPQ ) ) - AAPP = AAPP*DSQRT( ONE-T*AQOAP* - + AAPQ ) + MXSINJ = MAX( MXSINJ, DABS( SN ) ) + SVA( q ) = AAQQ*DSQRT( MAX( ZERO, + $ ONE+T*APOAQ*AAPQ ) ) + AAPP = AAPP*DSQRT( MAX( ZERO, + $ ONE-T*AQOAP*AAPQ ) ) * APOAQ = WORK( p ) / WORK( q ) AQOAP = WORK( q ) / WORK( p ) @@ -1177,26 +1274,26 @@ WORK( p ) = WORK( p )*CS WORK( q ) = WORK( q )*CS CALL DROTM( M, A( 1, p ), 1, - + A( 1, q ), 1, - + FASTR ) + $ A( 1, q ), 1, + $ FASTR ) IF( RSVEC )CALL DROTM( MVL, - + V( 1, p ), 1, V( 1, q ), - + 1, FASTR ) + $ V( 1, p ), 1, V( 1, q ), + $ 1, FASTR ) ELSE CALL DAXPY( M, -T*AQOAP, - + A( 1, q ), 1, - + A( 1, p ), 1 ) + $ A( 1, q ), 1, + $ A( 1, p ), 1 ) CALL DAXPY( M, CS*SN*APOAQ, - + A( 1, p ), 1, - + A( 1, q ), 1 ) + $ A( 1, p ), 1, + $ A( 1, q ), 1 ) IF( RSVEC ) THEN CALL DAXPY( MVL, -T*AQOAP, - + V( 1, q ), 1, - + V( 1, p ), 1 ) + $ V( 1, q ), 1, + $ V( 1, p ), 1 ) CALL DAXPY( MVL, - + CS*SN*APOAQ, - + V( 1, p ), 1, - + V( 1, q ), 1 ) + $ CS*SN*APOAQ, + $ V( 1, p ), 1, + $ V( 1, q ), 1 ) END IF WORK( p ) = WORK( p )*CS WORK( q ) = WORK( q ) / CS @@ -1204,61 +1301,61 @@ ELSE IF( WORK( q ).GE.ONE ) THEN CALL DAXPY( M, T*APOAQ, - + A( 1, p ), 1, - + A( 1, q ), 1 ) + $ A( 1, p ), 1, + $ A( 1, q ), 1 ) CALL DAXPY( M, -CS*SN*AQOAP, - + A( 1, q ), 1, - + A( 1, p ), 1 ) + $ A( 1, q ), 1, + $ A( 1, p ), 1 ) IF( RSVEC ) THEN CALL DAXPY( MVL, T*APOAQ, - + V( 1, p ), 1, - + V( 1, q ), 1 ) + $ V( 1, p ), 1, + $ V( 1, q ), 1 ) CALL DAXPY( MVL, - + -CS*SN*AQOAP, - + V( 1, q ), 1, - + V( 1, p ), 1 ) + $ -CS*SN*AQOAP, + $ V( 1, q ), 1, + $ V( 1, p ), 1 ) END IF WORK( p ) = WORK( p ) / CS WORK( q ) = WORK( q )*CS ELSE IF( WORK( p ).GE.WORK( q ) ) - + THEN + $ THEN CALL DAXPY( M, -T*AQOAP, - + A( 1, q ), 1, - + A( 1, p ), 1 ) + $ A( 1, q ), 1, + $ A( 1, p ), 1 ) CALL DAXPY( M, CS*SN*APOAQ, - + A( 1, p ), 1, - + A( 1, q ), 1 ) + $ A( 1, p ), 1, + $ A( 1, q ), 1 ) WORK( p ) = WORK( p )*CS WORK( q ) = WORK( q ) / CS IF( RSVEC ) THEN CALL DAXPY( MVL, - + -T*AQOAP, - + V( 1, q ), 1, - + V( 1, p ), 1 ) + $ -T*AQOAP, + $ V( 1, q ), 1, + $ V( 1, p ), 1 ) CALL DAXPY( MVL, - + CS*SN*APOAQ, - + V( 1, p ), 1, - + V( 1, q ), 1 ) + $ CS*SN*APOAQ, + $ V( 1, p ), 1, + $ V( 1, q ), 1 ) END IF ELSE CALL DAXPY( M, T*APOAQ, - + A( 1, p ), 1, - + A( 1, q ), 1 ) + $ A( 1, p ), 1, + $ A( 1, q ), 1 ) CALL DAXPY( M, - + -CS*SN*AQOAP, - + A( 1, q ), 1, - + A( 1, p ), 1 ) + $ -CS*SN*AQOAP, + $ A( 1, q ), 1, + $ A( 1, p ), 1 ) WORK( p ) = WORK( p ) / CS WORK( q ) = WORK( q )*CS IF( RSVEC ) THEN CALL DAXPY( MVL, - + T*APOAQ, V( 1, p ), - + 1, V( 1, q ), 1 ) + $ T*APOAQ, V( 1, p ), + $ 1, V( 1, q ), 1 ) CALL DAXPY( MVL, - + -CS*SN*AQOAP, - + V( 1, q ), 1, - + V( 1, p ), 1 ) + $ -CS*SN*AQOAP, + $ V( 1, q ), 1, + $ V( 1, p ), 1 ) END IF END IF END IF @@ -1268,40 +1365,40 @@ ELSE IF( AAPP.GT.AAQQ ) THEN CALL DCOPY( M, A( 1, p ), 1, - + WORK( N+1 ), 1 ) + $ WORK( N+1 ), 1 ) CALL DLASCL( 'G', 0, 0, AAPP, ONE, - + M, 1, WORK( N+1 ), LDA, - + IERR ) + $ M, 1, WORK( N+1 ), LDA, + $ IERR ) CALL DLASCL( 'G', 0, 0, AAQQ, ONE, - + M, 1, A( 1, q ), LDA, - + IERR ) + $ M, 1, A( 1, q ), LDA, + $ IERR ) TEMP1 = -AAPQ*WORK( p ) / WORK( q ) CALL DAXPY( M, TEMP1, WORK( N+1 ), - + 1, A( 1, q ), 1 ) + $ 1, A( 1, q ), 1 ) CALL DLASCL( 'G', 0, 0, ONE, AAQQ, - + M, 1, A( 1, q ), LDA, - + IERR ) - SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, - + ONE-AAPQ*AAPQ ) ) - MXSINJ = DMAX1( MXSINJ, SFMIN ) + $ M, 1, A( 1, q ), LDA, + $ IERR ) + SVA( q ) = AAQQ*DSQRT( MAX( ZERO, + $ ONE-AAPQ*AAPQ ) ) + MXSINJ = MAX( MXSINJ, SFMIN ) ELSE CALL DCOPY( M, A( 1, q ), 1, - + WORK( N+1 ), 1 ) + $ WORK( N+1 ), 1 ) CALL DLASCL( 'G', 0, 0, AAQQ, ONE, - + M, 1, WORK( N+1 ), LDA, - + IERR ) + $ M, 1, WORK( N+1 ), LDA, + $ IERR ) CALL DLASCL( 'G', 0, 0, AAPP, ONE, - + M, 1, A( 1, p ), LDA, - + IERR ) + $ M, 1, A( 1, p ), LDA, + $ IERR ) TEMP1 = -AAPQ*WORK( q ) / WORK( p ) CALL DAXPY( M, TEMP1, WORK( N+1 ), - + 1, A( 1, p ), 1 ) + $ 1, A( 1, p ), 1 ) CALL DLASCL( 'G', 0, 0, ONE, AAPP, - + M, 1, A( 1, p ), LDA, - + IERR ) - SVA( p ) = AAPP*DSQRT( DMAX1( ZERO, - + ONE-AAPQ*AAPQ ) ) - MXSINJ = DMAX1( MXSINJ, SFMIN ) + $ M, 1, A( 1, p ), LDA, + $ IERR ) + SVA( p ) = AAPP*DSQRT( MAX( ZERO, + $ ONE-AAPQ*AAPQ ) ) + MXSINJ = MAX( MXSINJ, SFMIN ) END IF END IF * END IF ROTOK THEN ... ELSE @@ -1309,29 +1406,29 @@ * In the case of cancellation in updating SVA(q) * .. recompute SVA(q) IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) - + THEN + $ THEN IF( ( AAQQ.LT.ROOTBIG ) .AND. - + ( AAQQ.GT.ROOTSFMIN ) ) THEN + $ ( AAQQ.GT.ROOTSFMIN ) ) THEN SVA( q ) = DNRM2( M, A( 1, q ), 1 )* - + WORK( q ) + $ WORK( q ) ELSE T = ZERO - AAQQ = ZERO + AAQQ = ONE CALL DLASSQ( M, A( 1, q ), 1, T, - + AAQQ ) + $ AAQQ ) SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) END IF END IF IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN IF( ( AAPP.LT.ROOTBIG ) .AND. - + ( AAPP.GT.ROOTSFMIN ) ) THEN + $ ( AAPP.GT.ROOTSFMIN ) ) THEN AAPP = DNRM2( M, A( 1, p ), 1 )* - + WORK( p ) + $ WORK( p ) ELSE T = ZERO - AAPP = ZERO + AAPP = ONE CALL DLASSQ( M, A( 1, p ), 1, T, - + AAPP ) + $ AAPP ) AAPP = T*DSQRT( AAPP )*WORK( p ) END IF SVA( p ) = AAPP @@ -1350,13 +1447,13 @@ END IF * IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) ) - + THEN + $ THEN SVA( p ) = AAPP NOTROT = 0 GO TO 2011 END IF IF( ( i.LE.SWBAND ) .AND. - + ( PSKIPPED.GT.ROWSKIP ) ) THEN + $ ( PSKIPPED.GT.ROWSKIP ) ) THEN AAPP = -AAPP NOTROT = 0 GO TO 2203 @@ -1371,7 +1468,7 @@ ELSE * IF( AAPP.EQ.ZERO )NOTROT = NOTROT + - + MIN0( jgl+KBL-1, N ) - jgl + 1 + $ MIN( jgl+KBL-1, N ) - jgl + 1 IF( AAPP.LT.ZERO )NOTROT = 0 * END IF @@ -1382,7 +1479,7 @@ * end of the jbc-loop 2011 CONTINUE *2011 bailed out of the jbc-loop - DO 2012 p = igl, MIN0( igl+KBL-1, N ) + DO 2012 p = igl, MIN( igl+KBL-1, N ) SVA( p ) = DABS( SVA( p ) ) 2012 CONTINUE *** @@ -1391,11 +1488,11 @@ * * .. update SVA(N) IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) ) - + THEN + $ THEN SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N ) ELSE T = ZERO - AAPP = ZERO + AAPP = ONE CALL DLASSQ( M, A( 1, N ), 1, T, AAPP ) SVA( N ) = T*DSQRT( AAPP )*WORK( N ) END IF @@ -1403,10 +1500,10 @@ * Additional steering devices * IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. - + ( ISWROT.LE.N ) ) )SWBAND = i + $ ( ISWROT.LE.N ) ) )SWBAND = i * IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )* - + TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN + $ TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN GO TO 1994 END IF * @@ -1446,12 +1543,12 @@ END IF IF( SVA( p ).NE.ZERO ) THEN N4 = N4 + 1 - IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1 + IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1 END IF 5991 CONTINUE IF( SVA( N ).NE.ZERO ) THEN N4 = N4 + 1 - IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1 + IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1 END IF * * Normalize the left singular vectors. @@ -1478,17 +1575,17 @@ END IF * * Undo scaling, if necessary (and possible). - IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / - + SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT. - + ( SFMIN / SCALE ) ) ) ) THEN + IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / SKL) ) ) + $ .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( MAX( N2, 1 ) ) .GT. + $ ( SFMIN / SKL) ) ) ) THEN DO 2400 p = 1, N - SVA( p ) = SCALE*SVA( p ) + SVA( P ) = SKL*SVA( P ) 2400 CONTINUE - SCALE = ONE + SKL= ONE END IF * - WORK( 1 ) = SCALE -* The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE + WORK( 1 ) = SKL +* The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE * then some of the singular values may overflow or underflow and * the spectrum is given in this factored representation. *