version 1.5, 2010/12/21 13:53:26
|
version 1.12, 2012/12/14 14:22:29
|
Line 1
|
Line 1
|
SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
*> \brief \b DGESVJ |
+ LDV, WORK, LWORK, INFO ) |
|
* |
* |
* -- LAPACK routine (version 3.3.0) -- |
* =========== DOCUMENTATION =========== |
* |
* |
* -- Contributed by Zlatko Drmac of the University of Zagreb and -- |
* Online html documentation available at |
* -- Kresimir Veselic of the Fernuniversitaet Hagen -- |
* http://www.netlib.org/lapack/explore-html/ |
* November 2010 |
|
* |
* |
|
*> \htmlonly |
|
*> Download DGESVJ + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f"> |
|
*> [TXT]</a> |
|
*> \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. |
|
*> \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 .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). |
|
*> \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 .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. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] MV |
|
*> \verbatim |
|
*> MV is 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. |
|
*> \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 .GE. 1. |
|
*> If JOBV .EQ. 'V', then LDV .GE. max(1,N). |
|
*> If JOBV .EQ. 'A', then LDV .GE. max(1,MV) . |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] WORK |
|
*> \verbatim |
|
*> WORK is 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 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 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. |
|
*> \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 September 2012 |
|
* |
|
*> \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.4.2) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
|
* September 2012 |
* |
* |
* 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 .. |
* .. Scalar Arguments .. |
INTEGER INFO, LDA, LDV, LWORK, M, MV, N |
INTEGER INFO, LDA, LDV, LWORK, M, MV, N |
CHARACTER*1 JOBA, JOBU, JOBV |
CHARACTER*1 JOBA, JOBU, JOBV |
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ), |
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 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 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 .. |
* .. Local Parameters .. |
DOUBLE PRECISION ZERO, HALF, ONE, TWO |
DOUBLE PRECISION ZERO, HALF, ONE |
PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, |
PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0) |
+ TWO = 2.0D0 ) |
|
INTEGER NSWEEP |
INTEGER NSWEEP |
PARAMETER ( NSWEEP = 30 ) |
PARAMETER ( NSWEEP = 30 ) |
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG, |
DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG, |
+ BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ, |
$ BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ, |
+ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, |
$ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, |
+ SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, |
$ SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, |
+ THSIGN, TOL |
$ THSIGN, TOL |
INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, |
INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, |
+ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, |
$ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, |
+ N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, |
$ N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, |
+ SWBAND |
$ SWBAND |
LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK, |
LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK, |
+ RSVEC, UCTOL, UPPER |
$ RSVEC, UCTOL, UPPER |
* .. |
* .. |
* .. Local Arrays .. |
* .. Local Arrays .. |
DOUBLE PRECISION FASTR( 5 ) |
DOUBLE PRECISION FASTR( 5 ) |
Line 327
|
Line 424
|
ELSE IF( MV.LT.0 ) THEN |
ELSE IF( MV.LT.0 ) THEN |
INFO = -9 |
INFO = -9 |
ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR. |
ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR. |
+ ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN |
$ ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN |
INFO = -11 |
INFO = -11 |
ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN |
ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN |
INFO = -12 |
INFO = -12 |
Line 383
|
Line 480
|
ROOTTOL = DSQRT( TOL ) |
ROOTTOL = DSQRT( TOL ) |
* |
* |
IF( DBLE( M )*EPSLN.GE.ONE ) THEN |
IF( DBLE( M )*EPSLN.GE.ONE ) THEN |
INFO = -5 |
INFO = -4 |
CALL XERBLA( 'DGESVJ', -INFO ) |
CALL XERBLA( 'DGESVJ', -INFO ) |
RETURN |
RETURN |
END IF |
END IF |
Line 518
|
Line 615
|
* |
* |
IF( N.EQ.1 ) THEN |
IF( N.EQ.1 ) THEN |
IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1, |
IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1, |
+ A( 1, 1 ), LDA, IERR ) |
$ A( 1, 1 ), LDA, IERR ) |
WORK( 1 ) = ONE / SKL |
WORK( 1 ) = ONE / SKL |
IF( SVA( 1 ).GE.SFMIN ) THEN |
IF( SVA( 1 ).GE.SFMIN ) THEN |
WORK( 2 ) = ONE |
WORK( 2 ) = ONE |
Line 538
|
Line 635
|
SN = DSQRT( SFMIN / EPSLN ) |
SN = DSQRT( SFMIN / EPSLN ) |
TEMP1 = DSQRT( BIG / DBLE( N ) ) |
TEMP1 = DSQRT( BIG / DBLE( N ) ) |
IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. |
IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. |
+ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN |
$ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN |
TEMP1 = DMIN1( BIG, TEMP1 / AAPP ) |
TEMP1 = DMIN1( BIG, TEMP1 / AAPP ) |
* AAQQ = AAQQ*TEMP1 |
* AAQQ = AAQQ*TEMP1 |
* AAPP = AAPP*TEMP1 |
* AAPP = AAPP*TEMP1 |
Line 638
|
Line 735
|
* [+ + x x] [x x]. [x x] |
* [+ + x x] [x x]. [x x] |
* |
* |
CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA, |
CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA, |
+ WORK( N34+1 ), SVA( N34+1 ), MVL, |
$ WORK( N34+1 ), SVA( N34+1 ), MVL, |
+ V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL, |
$ V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL, |
+ 2, WORK( N+1 ), LWORK-N, IERR ) |
$ 2, WORK( N+1 ), LWORK-N, IERR ) |
* |
* |
CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA, |
CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA, |
+ WORK( N2+1 ), SVA( N2+1 ), MVL, |
$ WORK( N2+1 ), SVA( N2+1 ), MVL, |
+ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2, |
$ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2, |
+ WORK( N+1 ), LWORK-N, IERR ) |
$ WORK( N+1 ), LWORK-N, IERR ) |
* |
* |
CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA, |
CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA, |
+ WORK( N2+1 ), SVA( N2+1 ), MVL, |
$ WORK( N2+1 ), SVA( N2+1 ), MVL, |
+ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, |
$ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, |
+ WORK( N+1 ), LWORK-N, IERR ) |
$ WORK( N+1 ), LWORK-N, IERR ) |
* |
* |
CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA, |
CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA, |
+ WORK( N4+1 ), SVA( N4+1 ), MVL, |
$ WORK( N4+1 ), SVA( N4+1 ), MVL, |
+ V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1, |
$ V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1, |
+ WORK( N+1 ), LWORK-N, IERR ) |
$ WORK( N+1 ), LWORK-N, IERR ) |
* |
* |
CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV, |
CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV, |
+ EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N, |
$ EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N, |
+ IERR ) |
$ IERR ) |
* |
* |
CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V, |
CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V, |
+ LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ), |
$ LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ), |
+ LWORK-N, IERR ) |
$ LWORK-N, IERR ) |
* |
* |
* |
* |
ELSE IF( UPPER ) THEN |
ELSE IF( UPPER ) THEN |
* |
* |
* |
* |
CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV, |
CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV, |
+ EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N, |
$ EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N, |
+ IERR ) |
$ IERR ) |
* |
* |
CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ), |
CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ), |
+ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV, |
$ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV, |
+ EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N, |
$ EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N, |
+ IERR ) |
$ IERR ) |
* |
* |
CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V, |
CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V, |
+ LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ), |
$ LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ), |
+ LWORK-N, IERR ) |
$ LWORK-N, IERR ) |
* |
* |
CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA, |
CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA, |
+ WORK( N2+1 ), SVA( N2+1 ), MVL, |
$ WORK( N2+1 ), SVA( N2+1 ), MVL, |
+ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, |
$ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, |
+ WORK( N+1 ), LWORK-N, IERR ) |
$ WORK( N+1 ), LWORK-N, IERR ) |
|
|
END IF |
END IF |
* |
* |
Line 725
|
Line 822
|
IF( p.NE.q ) THEN |
IF( p.NE.q ) THEN |
CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) |
CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) |
IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, |
IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, |
+ V( 1, q ), 1 ) |
$ V( 1, q ), 1 ) |
TEMP1 = SVA( p ) |
TEMP1 = SVA( p ) |
SVA( p ) = SVA( q ) |
SVA( p ) = SVA( q ) |
SVA( q ) = TEMP1 |
SVA( q ) = TEMP1 |
Line 749
|
Line 846
|
* below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)". |
* below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)". |
* |
* |
IF( ( SVA( p ).LT.ROOTBIG ) .AND. |
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 ) |
SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p ) |
ELSE |
ELSE |
TEMP1 = ZERO |
TEMP1 = ZERO |
Line 777
|
Line 874
|
ROTOK = ( SMALL*AAPP ).LE.AAQQ |
ROTOK = ( SMALL*AAPP ).LE.AAQQ |
IF( AAPP.LT.( BIG / AAQQ ) ) THEN |
IF( AAPP.LT.( BIG / AAQQ ) ) THEN |
AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, |
AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, |
+ q ), 1 )*WORK( p )*WORK( q ) / |
$ q ), 1 )*WORK( p )*WORK( q ) / |
+ AAQQ ) / AAPP |
$ AAQQ ) / AAPP |
ELSE |
ELSE |
CALL DCOPY( M, A( 1, p ), 1, |
CALL DCOPY( M, A( 1, p ), 1, |
+ WORK( N+1 ), 1 ) |
$ WORK( N+1 ), 1 ) |
CALL DLASCL( 'G', 0, 0, AAPP, |
CALL DLASCL( 'G', 0, 0, AAPP, |
+ WORK( p ), M, 1, |
$ WORK( p ), M, 1, |
+ WORK( N+1 ), LDA, IERR ) |
$ WORK( N+1 ), LDA, IERR ) |
AAPQ = DDOT( M, WORK( N+1 ), 1, |
AAPQ = DDOT( M, WORK( N+1 ), 1, |
+ A( 1, q ), 1 )*WORK( q ) / AAQQ |
$ A( 1, q ), 1 )*WORK( q ) / AAQQ |
END IF |
END IF |
ELSE |
ELSE |
ROTOK = AAPP.LE.( AAQQ / SMALL ) |
ROTOK = AAPP.LE.( AAQQ / SMALL ) |
IF( AAPP.GT.( SMALL / AAQQ ) ) THEN |
IF( AAPP.GT.( SMALL / AAQQ ) ) THEN |
AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, |
AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, |
+ q ), 1 )*WORK( p )*WORK( q ) / |
$ q ), 1 )*WORK( p )*WORK( q ) / |
+ AAQQ ) / AAPP |
$ AAQQ ) / AAPP |
ELSE |
ELSE |
CALL DCOPY( M, A( 1, q ), 1, |
CALL DCOPY( M, A( 1, q ), 1, |
+ WORK( N+1 ), 1 ) |
$ WORK( N+1 ), 1 ) |
CALL DLASCL( 'G', 0, 0, AAQQ, |
CALL DLASCL( 'G', 0, 0, AAQQ, |
+ WORK( q ), M, 1, |
$ WORK( q ), M, 1, |
+ WORK( N+1 ), LDA, IERR ) |
$ WORK( N+1 ), LDA, IERR ) |
AAPQ = DDOT( M, WORK( N+1 ), 1, |
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 |
END IF |
END IF |
* |
* |
Line 824
|
Line 921
|
* |
* |
AQOAP = AAQQ / AAPP |
AQOAP = AAQQ / AAPP |
APOAQ = AAPP / AAQQ |
APOAQ = AAPP / AAQQ |
THETA = -HALF*DABS( AQOAP-APOAQ ) / |
THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ |
+ AAPQ |
|
* |
* |
IF( DABS( THETA ).GT.BIGTHETA ) THEN |
IF( DABS( THETA ).GT.BIGTHETA ) THEN |
* |
* |
T = HALF / THETA |
T = HALF / THETA |
FASTR( 3 ) = T*WORK( p ) / WORK( q ) |
FASTR( 3 ) = T*WORK( p ) / WORK( q ) |
FASTR( 4 ) = -T*WORK( q ) / |
FASTR( 4 ) = -T*WORK( q ) / |
+ WORK( p ) |
$ WORK( p ) |
CALL DROTM( M, A( 1, p ), 1, |
CALL DROTM( M, A( 1, p ), 1, |
+ A( 1, q ), 1, FASTR ) |
$ A( 1, q ), 1, FASTR ) |
IF( RSVEC )CALL DROTM( MVL, |
IF( RSVEC )CALL DROTM( MVL, |
+ V( 1, p ), 1, |
$ V( 1, p ), 1, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ FASTR ) |
$ FASTR ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
+ ONE+T*APOAQ*AAPQ ) ) |
$ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
+ ONE-T*AQOAP*AAPQ ) ) |
$ ONE-T*AQOAP*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( T ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( T ) ) |
* |
* |
ELSE |
ELSE |
Line 851
|
Line 947
|
* |
* |
THSIGN = -DSIGN( ONE, AAPQ ) |
THSIGN = -DSIGN( ONE, AAPQ ) |
T = ONE / ( THETA+THSIGN* |
T = ONE / ( THETA+THSIGN* |
+ DSQRT( ONE+THETA*THETA ) ) |
$ DSQRT( ONE+THETA*THETA ) ) |
CS = DSQRT( ONE / ( ONE+T*T ) ) |
CS = DSQRT( ONE / ( ONE+T*T ) ) |
SN = T*CS |
SN = T*CS |
* |
* |
MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
+ ONE+T*APOAQ*AAPQ ) ) |
$ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
+ ONE-T*AQOAP*AAPQ ) ) |
$ ONE-T*AQOAP*AAPQ ) ) |
* |
* |
APOAQ = WORK( p ) / WORK( q ) |
APOAQ = WORK( p ) / WORK( q ) |
AQOAP = WORK( q ) / WORK( p ) |
AQOAP = WORK( q ) / WORK( p ) |
Line 870
|
Line 966
|
WORK( p ) = WORK( p )*CS |
WORK( p ) = WORK( p )*CS |
WORK( q ) = WORK( q )*CS |
WORK( q ) = WORK( q )*CS |
CALL DROTM( M, A( 1, p ), 1, |
CALL DROTM( M, A( 1, p ), 1, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ FASTR ) |
$ FASTR ) |
IF( RSVEC )CALL DROTM( MVL, |
IF( RSVEC )CALL DROTM( MVL, |
+ V( 1, p ), 1, V( 1, q ), |
$ V( 1, p ), 1, V( 1, q ), |
+ 1, FASTR ) |
$ 1, FASTR ) |
ELSE |
ELSE |
CALL DAXPY( M, -T*AQOAP, |
CALL DAXPY( M, -T*AQOAP, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ A( 1, p ), 1 ) |
$ A( 1, p ), 1 ) |
CALL DAXPY( M, CS*SN*APOAQ, |
CALL DAXPY( M, CS*SN*APOAQ, |
+ A( 1, p ), 1, |
$ A( 1, p ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
WORK( p ) = WORK( p )*CS |
WORK( p ) = WORK( p )*CS |
WORK( q ) = WORK( q ) / CS |
WORK( q ) = WORK( q ) / CS |
IF( RSVEC ) THEN |
IF( RSVEC ) THEN |
CALL DAXPY( MVL, -T*AQOAP, |
CALL DAXPY( MVL, -T*AQOAP, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ V( 1, p ), 1 ) |
$ V( 1, p ), 1 ) |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ CS*SN*APOAQ, |
$ CS*SN*APOAQ, |
+ V( 1, p ), 1, |
$ V( 1, p ), 1, |
+ V( 1, q ), 1 ) |
$ V( 1, q ), 1 ) |
END IF |
END IF |
END IF |
END IF |
ELSE |
ELSE |
IF( WORK( q ).GE.ONE ) THEN |
IF( WORK( q ).GE.ONE ) THEN |
CALL DAXPY( M, T*APOAQ, |
CALL DAXPY( M, T*APOAQ, |
+ A( 1, p ), 1, |
$ A( 1, p ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
CALL DAXPY( M, -CS*SN*AQOAP, |
CALL DAXPY( M, -CS*SN*AQOAP, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ A( 1, p ), 1 ) |
$ A( 1, p ), 1 ) |
WORK( p ) = WORK( p ) / CS |
WORK( p ) = WORK( p ) / CS |
WORK( q ) = WORK( q )*CS |
WORK( q ) = WORK( q )*CS |
IF( RSVEC ) THEN |
IF( RSVEC ) THEN |
CALL DAXPY( MVL, T*APOAQ, |
CALL DAXPY( MVL, T*APOAQ, |
+ V( 1, p ), 1, |
$ V( 1, p ), 1, |
+ V( 1, q ), 1 ) |
$ V( 1, q ), 1 ) |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ -CS*SN*AQOAP, |
$ -CS*SN*AQOAP, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ V( 1, p ), 1 ) |
$ V( 1, p ), 1 ) |
END IF |
END IF |
ELSE |
ELSE |
IF( WORK( p ).GE.WORK( q ) ) |
IF( WORK( p ).GE.WORK( q ) ) |
+ THEN |
$ THEN |
CALL DAXPY( M, -T*AQOAP, |
CALL DAXPY( M, -T*AQOAP, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ A( 1, p ), 1 ) |
$ A( 1, p ), 1 ) |
CALL DAXPY( M, CS*SN*APOAQ, |
CALL DAXPY( M, CS*SN*APOAQ, |
+ A( 1, p ), 1, |
$ A( 1, p ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
WORK( p ) = WORK( p )*CS |
WORK( p ) = WORK( p )*CS |
WORK( q ) = WORK( q ) / CS |
WORK( q ) = WORK( q ) / CS |
IF( RSVEC ) THEN |
IF( RSVEC ) THEN |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ -T*AQOAP, |
$ -T*AQOAP, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ V( 1, p ), 1 ) |
$ V( 1, p ), 1 ) |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ CS*SN*APOAQ, |
$ CS*SN*APOAQ, |
+ V( 1, p ), 1, |
$ V( 1, p ), 1, |
+ V( 1, q ), 1 ) |
$ V( 1, q ), 1 ) |
END IF |
END IF |
ELSE |
ELSE |
CALL DAXPY( M, T*APOAQ, |
CALL DAXPY( M, T*APOAQ, |
+ A( 1, p ), 1, |
$ A( 1, p ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
CALL DAXPY( M, |
CALL DAXPY( M, |
+ -CS*SN*AQOAP, |
$ -CS*SN*AQOAP, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ A( 1, p ), 1 ) |
$ A( 1, p ), 1 ) |
WORK( p ) = WORK( p ) / CS |
WORK( p ) = WORK( p ) / CS |
WORK( q ) = WORK( q )*CS |
WORK( q ) = WORK( q )*CS |
IF( RSVEC ) THEN |
IF( RSVEC ) THEN |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ T*APOAQ, V( 1, p ), |
$ T*APOAQ, V( 1, p ), |
+ 1, V( 1, q ), 1 ) |
$ 1, V( 1, q ), 1 ) |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ -CS*SN*AQOAP, |
$ -CS*SN*AQOAP, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ V( 1, p ), 1 ) |
$ V( 1, p ), 1 ) |
END IF |
END IF |
END IF |
END IF |
END IF |
END IF |
Line 961
|
Line 1057
|
ELSE |
ELSE |
* .. have to use modified Gram-Schmidt like transformation |
* .. have to use modified Gram-Schmidt like transformation |
CALL DCOPY( M, A( 1, p ), 1, |
CALL DCOPY( M, A( 1, p ), 1, |
+ WORK( N+1 ), 1 ) |
$ WORK( N+1 ), 1 ) |
CALL DLASCL( 'G', 0, 0, AAPP, ONE, M, |
CALL DLASCL( 'G', 0, 0, AAPP, ONE, M, |
+ 1, WORK( N+1 ), LDA, |
$ 1, WORK( N+1 ), LDA, |
+ IERR ) |
$ IERR ) |
CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M, |
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 ) |
TEMP1 = -AAPQ*WORK( p ) / WORK( q ) |
CALL DAXPY( M, TEMP1, WORK( N+1 ), 1, |
CALL DAXPY( M, TEMP1, WORK( N+1 ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M, |
CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M, |
+ 1, A( 1, q ), LDA, IERR ) |
$ 1, A( 1, q ), LDA, IERR ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
+ ONE-AAPQ*AAPQ ) ) |
$ ONE-AAPQ*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
END IF |
END IF |
* END IF ROTOK THEN ... ELSE |
* END IF ROTOK THEN ... ELSE |
Line 982
|
Line 1078
|
* recompute SVA(q), SVA(p). |
* recompute SVA(q), SVA(p). |
* |
* |
IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) |
IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) |
+ THEN |
$ THEN |
IF( ( AAQQ.LT.ROOTBIG ) .AND. |
IF( ( AAQQ.LT.ROOTBIG ) .AND. |
+ ( AAQQ.GT.ROOTSFMIN ) ) THEN |
$ ( AAQQ.GT.ROOTSFMIN ) ) THEN |
SVA( q ) = DNRM2( M, A( 1, q ), 1 )* |
SVA( q ) = DNRM2( M, A( 1, q ), 1 )* |
+ WORK( q ) |
$ WORK( q ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAQQ = ONE |
AAQQ = ONE |
CALL DLASSQ( M, A( 1, q ), 1, T, |
CALL DLASSQ( M, A( 1, q ), 1, T, |
+ AAQQ ) |
$ AAQQ ) |
SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) |
SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) |
END IF |
END IF |
END IF |
END IF |
IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN |
IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN |
IF( ( AAPP.LT.ROOTBIG ) .AND. |
IF( ( AAPP.LT.ROOTBIG ) .AND. |
+ ( AAPP.GT.ROOTSFMIN ) ) THEN |
$ ( AAPP.GT.ROOTSFMIN ) ) THEN |
AAPP = DNRM2( M, A( 1, p ), 1 )* |
AAPP = DNRM2( M, A( 1, p ), 1 )* |
+ WORK( p ) |
$ WORK( p ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAPP = ONE |
AAPP = ONE |
CALL DLASSQ( M, A( 1, p ), 1, T, |
CALL DLASSQ( M, A( 1, p ), 1, T, |
+ AAPP ) |
$ AAPP ) |
AAPP = T*DSQRT( AAPP )*WORK( p ) |
AAPP = T*DSQRT( AAPP )*WORK( p ) |
END IF |
END IF |
SVA( p ) = AAPP |
SVA( p ) = AAPP |
Line 1023
|
Line 1119
|
END IF |
END IF |
* |
* |
IF( ( i.LE.SWBAND ) .AND. |
IF( ( i.LE.SWBAND ) .AND. |
+ ( PSKIPPED.GT.ROWSKIP ) ) THEN |
$ ( PSKIPPED.GT.ROWSKIP ) ) THEN |
IF( ir1.EQ.0 )AAPP = -AAPP |
IF( ir1.EQ.0 )AAPP = -AAPP |
NOTROT = 0 |
NOTROT = 0 |
GO TO 2103 |
GO TO 2103 |
Line 1040
|
Line 1136
|
ELSE |
ELSE |
SVA( p ) = AAPP |
SVA( p ) = AAPP |
IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) ) |
IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) ) |
+ NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p |
$ NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p |
END IF |
END IF |
* |
* |
2001 CONTINUE |
2001 CONTINUE |
Line 1085
|
Line 1181
|
END IF |
END IF |
IF( AAPP.LT.( BIG / AAQQ ) ) THEN |
IF( AAPP.LT.( BIG / AAQQ ) ) THEN |
AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, |
AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, |
+ q ), 1 )*WORK( p )*WORK( q ) / |
$ q ), 1 )*WORK( p )*WORK( q ) / |
+ AAQQ ) / AAPP |
$ AAQQ ) / AAPP |
ELSE |
ELSE |
CALL DCOPY( M, A( 1, p ), 1, |
CALL DCOPY( M, A( 1, p ), 1, |
+ WORK( N+1 ), 1 ) |
$ WORK( N+1 ), 1 ) |
CALL DLASCL( 'G', 0, 0, AAPP, |
CALL DLASCL( 'G', 0, 0, AAPP, |
+ WORK( p ), M, 1, |
$ WORK( p ), M, 1, |
+ WORK( N+1 ), LDA, IERR ) |
$ WORK( N+1 ), LDA, IERR ) |
AAPQ = DDOT( M, WORK( N+1 ), 1, |
AAPQ = DDOT( M, WORK( N+1 ), 1, |
+ A( 1, q ), 1 )*WORK( q ) / AAQQ |
$ A( 1, q ), 1 )*WORK( q ) / AAQQ |
END IF |
END IF |
ELSE |
ELSE |
IF( AAPP.GE.AAQQ ) THEN |
IF( AAPP.GE.AAQQ ) THEN |
Line 1104
|
Line 1200
|
END IF |
END IF |
IF( AAPP.GT.( SMALL / AAQQ ) ) THEN |
IF( AAPP.GT.( SMALL / AAQQ ) ) THEN |
AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, |
AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1, |
+ q ), 1 )*WORK( p )*WORK( q ) / |
$ q ), 1 )*WORK( p )*WORK( q ) / |
+ AAQQ ) / AAPP |
$ AAQQ ) / AAPP |
ELSE |
ELSE |
CALL DCOPY( M, A( 1, q ), 1, |
CALL DCOPY( M, A( 1, q ), 1, |
+ WORK( N+1 ), 1 ) |
$ WORK( N+1 ), 1 ) |
CALL DLASCL( 'G', 0, 0, AAQQ, |
CALL DLASCL( 'G', 0, 0, AAQQ, |
+ WORK( q ), M, 1, |
$ WORK( q ), M, 1, |
+ WORK( N+1 ), LDA, IERR ) |
$ WORK( N+1 ), LDA, IERR ) |
AAPQ = DDOT( M, WORK( N+1 ), 1, |
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 |
END IF |
END IF |
* |
* |
Line 1131
|
Line 1227
|
* |
* |
AQOAP = AAQQ / AAPP |
AQOAP = AAQQ / AAPP |
APOAQ = AAPP / AAQQ |
APOAQ = AAPP / AAQQ |
THETA = -HALF*DABS( AQOAP-APOAQ ) / |
THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ |
+ AAPQ |
|
IF( AAQQ.GT.AAPP0 )THETA = -THETA |
IF( AAQQ.GT.AAPP0 )THETA = -THETA |
* |
* |
IF( DABS( THETA ).GT.BIGTHETA ) THEN |
IF( DABS( THETA ).GT.BIGTHETA ) THEN |
T = HALF / THETA |
T = HALF / THETA |
FASTR( 3 ) = T*WORK( p ) / WORK( q ) |
FASTR( 3 ) = T*WORK( p ) / WORK( q ) |
FASTR( 4 ) = -T*WORK( q ) / |
FASTR( 4 ) = -T*WORK( q ) / |
+ WORK( p ) |
$ WORK( p ) |
CALL DROTM( M, A( 1, p ), 1, |
CALL DROTM( M, A( 1, p ), 1, |
+ A( 1, q ), 1, FASTR ) |
$ A( 1, q ), 1, FASTR ) |
IF( RSVEC )CALL DROTM( MVL, |
IF( RSVEC )CALL DROTM( MVL, |
+ V( 1, p ), 1, |
$ V( 1, p ), 1, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ FASTR ) |
$ FASTR ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
+ ONE+T*APOAQ*AAPQ ) ) |
$ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
+ ONE-T*AQOAP*AAPQ ) ) |
$ ONE-T*AQOAP*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( T ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( T ) ) |
ELSE |
ELSE |
* |
* |
Line 1158
|
Line 1253
|
THSIGN = -DSIGN( ONE, AAPQ ) |
THSIGN = -DSIGN( ONE, AAPQ ) |
IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN |
IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN |
T = ONE / ( THETA+THSIGN* |
T = ONE / ( THETA+THSIGN* |
+ DSQRT( ONE+THETA*THETA ) ) |
$ DSQRT( ONE+THETA*THETA ) ) |
CS = DSQRT( ONE / ( ONE+T*T ) ) |
CS = DSQRT( ONE / ( ONE+T*T ) ) |
SN = T*CS |
SN = T*CS |
MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
+ ONE+T*APOAQ*AAPQ ) ) |
$ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
+ ONE-T*AQOAP*AAPQ ) ) |
$ ONE-T*AQOAP*AAPQ ) ) |
* |
* |
APOAQ = WORK( p ) / WORK( q ) |
APOAQ = WORK( p ) / WORK( q ) |
AQOAP = WORK( q ) / WORK( p ) |
AQOAP = WORK( q ) / WORK( p ) |
Line 1177
|
Line 1272
|
WORK( p ) = WORK( p )*CS |
WORK( p ) = WORK( p )*CS |
WORK( q ) = WORK( q )*CS |
WORK( q ) = WORK( q )*CS |
CALL DROTM( M, A( 1, p ), 1, |
CALL DROTM( M, A( 1, p ), 1, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ FASTR ) |
$ FASTR ) |
IF( RSVEC )CALL DROTM( MVL, |
IF( RSVEC )CALL DROTM( MVL, |
+ V( 1, p ), 1, V( 1, q ), |
$ V( 1, p ), 1, V( 1, q ), |
+ 1, FASTR ) |
$ 1, FASTR ) |
ELSE |
ELSE |
CALL DAXPY( M, -T*AQOAP, |
CALL DAXPY( M, -T*AQOAP, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ A( 1, p ), 1 ) |
$ A( 1, p ), 1 ) |
CALL DAXPY( M, CS*SN*APOAQ, |
CALL DAXPY( M, CS*SN*APOAQ, |
+ A( 1, p ), 1, |
$ A( 1, p ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
IF( RSVEC ) THEN |
IF( RSVEC ) THEN |
CALL DAXPY( MVL, -T*AQOAP, |
CALL DAXPY( MVL, -T*AQOAP, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ V( 1, p ), 1 ) |
$ V( 1, p ), 1 ) |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ CS*SN*APOAQ, |
$ CS*SN*APOAQ, |
+ V( 1, p ), 1, |
$ V( 1, p ), 1, |
+ V( 1, q ), 1 ) |
$ V( 1, q ), 1 ) |
END IF |
END IF |
WORK( p ) = WORK( p )*CS |
WORK( p ) = WORK( p )*CS |
WORK( q ) = WORK( q ) / CS |
WORK( q ) = WORK( q ) / CS |
Line 1204
|
Line 1299
|
ELSE |
ELSE |
IF( WORK( q ).GE.ONE ) THEN |
IF( WORK( q ).GE.ONE ) THEN |
CALL DAXPY( M, T*APOAQ, |
CALL DAXPY( M, T*APOAQ, |
+ A( 1, p ), 1, |
$ A( 1, p ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
CALL DAXPY( M, -CS*SN*AQOAP, |
CALL DAXPY( M, -CS*SN*AQOAP, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ A( 1, p ), 1 ) |
$ A( 1, p ), 1 ) |
IF( RSVEC ) THEN |
IF( RSVEC ) THEN |
CALL DAXPY( MVL, T*APOAQ, |
CALL DAXPY( MVL, T*APOAQ, |
+ V( 1, p ), 1, |
$ V( 1, p ), 1, |
+ V( 1, q ), 1 ) |
$ V( 1, q ), 1 ) |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ -CS*SN*AQOAP, |
$ -CS*SN*AQOAP, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ V( 1, p ), 1 ) |
$ V( 1, p ), 1 ) |
END IF |
END IF |
WORK( p ) = WORK( p ) / CS |
WORK( p ) = WORK( p ) / CS |
WORK( q ) = WORK( q )*CS |
WORK( q ) = WORK( q )*CS |
ELSE |
ELSE |
IF( WORK( p ).GE.WORK( q ) ) |
IF( WORK( p ).GE.WORK( q ) ) |
+ THEN |
$ THEN |
CALL DAXPY( M, -T*AQOAP, |
CALL DAXPY( M, -T*AQOAP, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ A( 1, p ), 1 ) |
$ A( 1, p ), 1 ) |
CALL DAXPY( M, CS*SN*APOAQ, |
CALL DAXPY( M, CS*SN*APOAQ, |
+ A( 1, p ), 1, |
$ A( 1, p ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
WORK( p ) = WORK( p )*CS |
WORK( p ) = WORK( p )*CS |
WORK( q ) = WORK( q ) / CS |
WORK( q ) = WORK( q ) / CS |
IF( RSVEC ) THEN |
IF( RSVEC ) THEN |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ -T*AQOAP, |
$ -T*AQOAP, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ V( 1, p ), 1 ) |
$ V( 1, p ), 1 ) |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ CS*SN*APOAQ, |
$ CS*SN*APOAQ, |
+ V( 1, p ), 1, |
$ V( 1, p ), 1, |
+ V( 1, q ), 1 ) |
$ V( 1, q ), 1 ) |
END IF |
END IF |
ELSE |
ELSE |
CALL DAXPY( M, T*APOAQ, |
CALL DAXPY( M, T*APOAQ, |
+ A( 1, p ), 1, |
$ A( 1, p ), 1, |
+ A( 1, q ), 1 ) |
$ A( 1, q ), 1 ) |
CALL DAXPY( M, |
CALL DAXPY( M, |
+ -CS*SN*AQOAP, |
$ -CS*SN*AQOAP, |
+ A( 1, q ), 1, |
$ A( 1, q ), 1, |
+ A( 1, p ), 1 ) |
$ A( 1, p ), 1 ) |
WORK( p ) = WORK( p ) / CS |
WORK( p ) = WORK( p ) / CS |
WORK( q ) = WORK( q )*CS |
WORK( q ) = WORK( q )*CS |
IF( RSVEC ) THEN |
IF( RSVEC ) THEN |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ T*APOAQ, V( 1, p ), |
$ T*APOAQ, V( 1, p ), |
+ 1, V( 1, q ), 1 ) |
$ 1, V( 1, q ), 1 ) |
CALL DAXPY( MVL, |
CALL DAXPY( MVL, |
+ -CS*SN*AQOAP, |
$ -CS*SN*AQOAP, |
+ V( 1, q ), 1, |
$ V( 1, q ), 1, |
+ V( 1, p ), 1 ) |
$ V( 1, p ), 1 ) |
END IF |
END IF |
END IF |
END IF |
END IF |
END IF |
Line 1268
|
Line 1363
|
ELSE |
ELSE |
IF( AAPP.GT.AAQQ ) THEN |
IF( AAPP.GT.AAQQ ) THEN |
CALL DCOPY( M, A( 1, p ), 1, |
CALL DCOPY( M, A( 1, p ), 1, |
+ WORK( N+1 ), 1 ) |
$ WORK( N+1 ), 1 ) |
CALL DLASCL( 'G', 0, 0, AAPP, ONE, |
CALL DLASCL( 'G', 0, 0, AAPP, ONE, |
+ M, 1, WORK( N+1 ), LDA, |
$ M, 1, WORK( N+1 ), LDA, |
+ IERR ) |
$ IERR ) |
CALL DLASCL( 'G', 0, 0, AAQQ, ONE, |
CALL DLASCL( 'G', 0, 0, AAQQ, ONE, |
+ M, 1, A( 1, q ), LDA, |
$ M, 1, A( 1, q ), LDA, |
+ IERR ) |
$ IERR ) |
TEMP1 = -AAPQ*WORK( p ) / WORK( q ) |
TEMP1 = -AAPQ*WORK( p ) / WORK( q ) |
CALL DAXPY( M, TEMP1, WORK( N+1 ), |
CALL DAXPY( M, TEMP1, WORK( N+1 ), |
+ 1, A( 1, q ), 1 ) |
$ 1, A( 1, q ), 1 ) |
CALL DLASCL( 'G', 0, 0, ONE, AAQQ, |
CALL DLASCL( 'G', 0, 0, ONE, AAQQ, |
+ M, 1, A( 1, q ), LDA, |
$ M, 1, A( 1, q ), LDA, |
+ IERR ) |
$ IERR ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
+ ONE-AAPQ*AAPQ ) ) |
$ ONE-AAPQ*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
ELSE |
ELSE |
CALL DCOPY( M, A( 1, q ), 1, |
CALL DCOPY( M, A( 1, q ), 1, |
+ WORK( N+1 ), 1 ) |
$ WORK( N+1 ), 1 ) |
CALL DLASCL( 'G', 0, 0, AAQQ, ONE, |
CALL DLASCL( 'G', 0, 0, AAQQ, ONE, |
+ M, 1, WORK( N+1 ), LDA, |
$ M, 1, WORK( N+1 ), LDA, |
+ IERR ) |
$ IERR ) |
CALL DLASCL( 'G', 0, 0, AAPP, ONE, |
CALL DLASCL( 'G', 0, 0, AAPP, ONE, |
+ M, 1, A( 1, p ), LDA, |
$ M, 1, A( 1, p ), LDA, |
+ IERR ) |
$ IERR ) |
TEMP1 = -AAPQ*WORK( q ) / WORK( p ) |
TEMP1 = -AAPQ*WORK( q ) / WORK( p ) |
CALL DAXPY( M, TEMP1, WORK( N+1 ), |
CALL DAXPY( M, TEMP1, WORK( N+1 ), |
+ 1, A( 1, p ), 1 ) |
$ 1, A( 1, p ), 1 ) |
CALL DLASCL( 'G', 0, 0, ONE, AAPP, |
CALL DLASCL( 'G', 0, 0, ONE, AAPP, |
+ M, 1, A( 1, p ), LDA, |
$ M, 1, A( 1, p ), LDA, |
+ IERR ) |
$ IERR ) |
SVA( p ) = AAPP*DSQRT( DMAX1( ZERO, |
SVA( p ) = AAPP*DSQRT( DMAX1( ZERO, |
+ ONE-AAPQ*AAPQ ) ) |
$ ONE-AAPQ*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
END IF |
END IF |
END IF |
END IF |
Line 1309
|
Line 1404
|
* In the case of cancellation in updating SVA(q) |
* In the case of cancellation in updating SVA(q) |
* .. recompute SVA(q) |
* .. recompute SVA(q) |
IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) |
IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) |
+ THEN |
$ THEN |
IF( ( AAQQ.LT.ROOTBIG ) .AND. |
IF( ( AAQQ.LT.ROOTBIG ) .AND. |
+ ( AAQQ.GT.ROOTSFMIN ) ) THEN |
$ ( AAQQ.GT.ROOTSFMIN ) ) THEN |
SVA( q ) = DNRM2( M, A( 1, q ), 1 )* |
SVA( q ) = DNRM2( M, A( 1, q ), 1 )* |
+ WORK( q ) |
$ WORK( q ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAQQ = ONE |
AAQQ = ONE |
CALL DLASSQ( M, A( 1, q ), 1, T, |
CALL DLASSQ( M, A( 1, q ), 1, T, |
+ AAQQ ) |
$ AAQQ ) |
SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) |
SVA( q ) = T*DSQRT( AAQQ )*WORK( q ) |
END IF |
END IF |
END IF |
END IF |
IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN |
IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN |
IF( ( AAPP.LT.ROOTBIG ) .AND. |
IF( ( AAPP.LT.ROOTBIG ) .AND. |
+ ( AAPP.GT.ROOTSFMIN ) ) THEN |
$ ( AAPP.GT.ROOTSFMIN ) ) THEN |
AAPP = DNRM2( M, A( 1, p ), 1 )* |
AAPP = DNRM2( M, A( 1, p ), 1 )* |
+ WORK( p ) |
$ WORK( p ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
AAPP = ONE |
AAPP = ONE |
CALL DLASSQ( M, A( 1, p ), 1, T, |
CALL DLASSQ( M, A( 1, p ), 1, T, |
+ AAPP ) |
$ AAPP ) |
AAPP = T*DSQRT( AAPP )*WORK( p ) |
AAPP = T*DSQRT( AAPP )*WORK( p ) |
END IF |
END IF |
SVA( p ) = AAPP |
SVA( p ) = AAPP |
Line 1350
|
Line 1445
|
END IF |
END IF |
* |
* |
IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) ) |
IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) ) |
+ THEN |
$ THEN |
SVA( p ) = AAPP |
SVA( p ) = AAPP |
NOTROT = 0 |
NOTROT = 0 |
GO TO 2011 |
GO TO 2011 |
END IF |
END IF |
IF( ( i.LE.SWBAND ) .AND. |
IF( ( i.LE.SWBAND ) .AND. |
+ ( PSKIPPED.GT.ROWSKIP ) ) THEN |
$ ( PSKIPPED.GT.ROWSKIP ) ) THEN |
AAPP = -AAPP |
AAPP = -AAPP |
NOTROT = 0 |
NOTROT = 0 |
GO TO 2203 |
GO TO 2203 |
Line 1371
|
Line 1466
|
ELSE |
ELSE |
* |
* |
IF( AAPP.EQ.ZERO )NOTROT = NOTROT + |
IF( AAPP.EQ.ZERO )NOTROT = NOTROT + |
+ MIN0( jgl+KBL-1, N ) - jgl + 1 |
$ MIN0( jgl+KBL-1, N ) - jgl + 1 |
IF( AAPP.LT.ZERO )NOTROT = 0 |
IF( AAPP.LT.ZERO )NOTROT = 0 |
* |
* |
END IF |
END IF |
Line 1391
|
Line 1486
|
* |
* |
* .. update SVA(N) |
* .. update SVA(N) |
IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) ) |
IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) ) |
+ THEN |
$ THEN |
SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N ) |
SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N ) |
ELSE |
ELSE |
T = ZERO |
T = ZERO |
Line 1403
|
Line 1498
|
* Additional steering devices |
* Additional steering devices |
* |
* |
IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. |
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 ) )* |
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 |
GO TO 1994 |
END IF |
END IF |
* |
* |
Line 1478
|
Line 1573
|
END IF |
END IF |
* |
* |
* Undo scaling, if necessary (and possible). |
* Undo scaling, if necessary (and possible). |
IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / |
IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / SKL) ) ) |
+ SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT. |
$ .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( MAX( N2, 1 ) ) .GT. |
+ ( SFMIN / SKL) ) ) ) THEN |
$ ( SFMIN / SKL) ) ) ) THEN |
DO 2400 p = 1, N |
DO 2400 p = 1, N |
SVA( p ) = SKL*SVA( p ) |
SVA( P ) = SKL*SVA( P ) |
2400 CONTINUE |
2400 CONTINUE |
SKL= ONE |
SKL= ONE |
END IF |
END IF |