version 1.12, 2012/12/14 14:22:29
|
version 1.21, 2023/08/07 08:38:50
|
Line 2
|
Line 2
|
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download DGESVJ + dependencies |
*> Download DGESVJ + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
* SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
* LDV, WORK, LWORK, INFO ) |
* LDV, WORK, LWORK, INFO ) |
* |
* |
* .. 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 |
Line 29
|
Line 29
|
* DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ), |
* DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ), |
* $ WORK( LWORK ) |
* $ WORK( LWORK ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 45
|
Line 45
|
*> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements |
*> 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 |
*> 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. |
*> 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 |
*> \endverbatim |
* |
* |
* Arguments: |
* Arguments: |
Line 52
|
Line 54
|
* |
* |
*> \param[in] JOBA |
*> \param[in] JOBA |
*> \verbatim |
*> \verbatim |
*> JOBA is CHARACTER* 1 |
*> JOBA is CHARACTER*1 |
*> Specifies the structure of A. |
*> Specifies the structure of A. |
*> = 'L': The input matrix A is lower triangular; |
*> = 'L': The input matrix A is lower triangular; |
*> = 'U': The input matrix A is upper triangular; |
*> = 'U': The input matrix A is upper triangular; |
Line 88
|
Line 90
|
*> JOBV is CHARACTER*1 |
*> JOBV is CHARACTER*1 |
*> Specifies whether to compute the right singular vectors, that |
*> Specifies whether to compute the right singular vectors, that |
*> is, the matrix V: |
*> is, the matrix V: |
*> = 'V' : the matrix V is computed and returned in the array V |
*> = 'V': the matrix V is computed and returned in the array V |
*> = 'A' : the Jacobi rotations are applied to the MV-by-N |
*> = 'A': the Jacobi rotations are applied to the MV-by-N |
*> array V. In other words, the right singular vector |
*> array V. In other words, the right singular vector |
*> matrix V is not computed explicitly, instead it is |
*> matrix V is not computed explicitly, instead it is |
*> applied to an MV-by-N matrix initially stored in the |
*> applied to an MV-by-N matrix initially stored in the |
*> first MV rows of V. |
*> first MV rows of V. |
*> = 'N' : the matrix V is not computed and the array V is not |
*> = 'N': the matrix V is not computed and the array V is not |
*> referenced |
*> referenced |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] M |
*> \param[in] M |
*> \verbatim |
*> \verbatim |
*> M is INTEGER |
*> M is INTEGER |
*> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0. |
*> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] N |
*> \param[in] N |
Line 116
|
Line 118
|
*> A is DOUBLE PRECISION array, dimension (LDA,N) |
*> A is DOUBLE PRECISION array, dimension (LDA,N) |
*> On entry, the M-by-N matrix A. |
*> On entry, the M-by-N matrix A. |
*> On exit : |
*> On exit : |
*> If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' : |
*> If JOBU = 'U' .OR. JOBU = 'C' : |
*> If INFO .EQ. 0 : |
*> If INFO = 0 : |
*> RANKA orthonormal columns of U are returned in the |
*> RANKA orthonormal columns of U are returned in the |
*> leading RANKA columns of the array A. Here RANKA <= N |
*> leading RANKA columns of the array A. Here RANKA <= N |
*> is the number of computed singular values of A that are |
*> is the number of computed singular values of A that are |
Line 127
|
Line 129
|
*> in the array WORK as RANKA=NINT(WORK(2)). Also see the |
*> in the array WORK as RANKA=NINT(WORK(2)). Also see the |
*> descriptions of SVA and WORK. The computed columns of U |
*> descriptions of SVA and WORK. The computed columns of U |
*> are mutually numerically orthogonal up to approximately |
*> are mutually numerically orthogonal up to approximately |
*> TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'), |
*> TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), |
*> see the description of JOBU. |
*> see the description of JOBU. |
*> If INFO .GT. 0 : |
*> If INFO > 0 : |
*> the procedure DGESVJ did not converge in the given number |
*> the procedure DGESVJ did not converge in the given number |
*> of iterations (sweeps). In that case, the computed |
*> of iterations (sweeps). In that case, the computed |
*> columns of U may not be orthogonal up to TOL. The output |
*> columns of U may not be orthogonal up to TOL. The output |
Line 138
|
Line 140
|
*> input matrix A in the sense that the residual |
*> input matrix A in the sense that the residual |
*> ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small. |
*> ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small. |
*> |
*> |
*> If JOBU .EQ. 'N' : |
*> If JOBU = 'N' : |
*> If INFO .EQ. 0 : |
*> If INFO = 0 : |
*> Note that the left singular vectors are 'for free' in the |
*> Note that the left singular vectors are 'for free' in the |
*> one-sided Jacobi SVD algorithm. However, if only the |
*> one-sided Jacobi SVD algorithm. However, if only the |
*> singular values are needed, the level of numerical |
*> singular values are needed, the level of numerical |
Line 148
|
Line 150
|
*> numerically orthogonal up to approximately M*EPS. Thus, |
*> numerically orthogonal up to approximately M*EPS. Thus, |
*> on exit, A contains the columns of U scaled with the |
*> on exit, A contains the columns of U scaled with the |
*> corresponding singular values. |
*> corresponding singular values. |
*> If INFO .GT. 0 : |
*> If INFO > 0 : |
*> the procedure DGESVJ did not converge in the given number |
*> the procedure DGESVJ did not converge in the given number |
*> of iterations (sweeps). |
*> of iterations (sweeps). |
*> \endverbatim |
*> \endverbatim |
Line 163
|
Line 165
|
*> \verbatim |
*> \verbatim |
*> SVA is DOUBLE PRECISION array, dimension (N) |
*> SVA is DOUBLE PRECISION array, dimension (N) |
*> On exit : |
*> On exit : |
*> If INFO .EQ. 0 : |
*> If INFO = 0 : |
*> depending on the value SCALE = WORK(1), we have: |
*> depending on the value SCALE = WORK(1), we have: |
*> If SCALE .EQ. ONE : |
*> If SCALE = ONE : |
*> SVA(1:N) contains the computed singular values of A. |
*> SVA(1:N) contains the computed singular values of A. |
*> During the computation SVA contains the Euclidean column |
*> During the computation SVA contains the Euclidean column |
*> norms of the iterated matrices in the array A. |
*> norms of the iterated matrices in the array A. |
Line 173
|
Line 175
|
*> The singular values of A are SCALE*SVA(1:N), and this |
*> The singular values of A are SCALE*SVA(1:N), and this |
*> factored representation is due to the fact that some of the |
*> factored representation is due to the fact that some of the |
*> singular values of A might underflow or overflow. |
*> singular values of A might underflow or overflow. |
*> If INFO .GT. 0 : |
*> If INFO > 0 : |
*> the procedure DGESVJ did not converge in the given number of |
*> the procedure DGESVJ did not converge in the given number of |
*> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate. |
*> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate. |
*> \endverbatim |
*> \endverbatim |
Line 181
|
Line 183
|
*> \param[in] MV |
*> \param[in] MV |
*> \verbatim |
*> \verbatim |
*> MV is INTEGER |
*> MV is INTEGER |
*> If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ |
*> 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. |
*> is applied to the first MV rows of V. See the description of JOBV. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 199
|
Line 201
|
*> \param[in] LDV |
*> \param[in] LDV |
*> \verbatim |
*> \verbatim |
*> LDV is INTEGER |
*> LDV is INTEGER |
*> The leading dimension of the array V, LDV .GE. 1. |
*> The leading dimension of the array V, LDV >= 1. |
*> If JOBV .EQ. 'V', then LDV .GE. max(1,N). |
*> If JOBV = 'V', then LDV >= max(1,N). |
*> If JOBV .EQ. 'A', then LDV .GE. max(1,MV) . |
*> If JOBV = 'A', then LDV >= max(1,MV) . |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in,out] WORK |
*> \param[in,out] WORK |
*> \verbatim |
*> \verbatim |
*> WORK is DOUBLE PRECISION array, dimension max(4,M+N). |
*> WORK is DOUBLE PRECISION array, dimension (LWORK) |
*> On entry : |
*> On entry : |
*> If JOBU .EQ. 'C' : |
*> If JOBU = 'C' : |
*> WORK(1) = CTOL, where CTOL defines the threshold for convergence. |
*> WORK(1) = CTOL, where CTOL defines the threshold for convergence. |
*> The process stops if all columns of A are mutually |
*> The process stops if all columns of A are mutually |
*> orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). |
*> orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). |
Line 228
|
Line 230
|
*> WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. |
*> WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. |
*> This is useful information in cases when DGESVJ did |
*> This is useful information in cases when DGESVJ did |
*> not converge, as it can be used to estimate whether |
*> not converge, as it can be used to estimate whether |
*> the output is stil useful and for post festum analysis. |
*> the output is still useful and for post festum analysis. |
*> WORK(6) = the largest absolute value over all sines of the |
*> WORK(6) = the largest absolute value over all sines of the |
*> Jacobi rotation angles in the last sweep. It can be |
*> Jacobi rotation angles in the last sweep. It can be |
*> useful for a post festum analysis. |
*> useful for a post festum analysis. |
Line 243
|
Line 245
|
*> \param[out] INFO |
*> \param[out] INFO |
*> \verbatim |
*> \verbatim |
*> INFO is INTEGER |
*> INFO is INTEGER |
*> = 0 : successful exit. |
*> = 0: successful exit. |
*> < 0 : if INFO = -i, then the i-th argument had an illegal value |
*> < 0: if INFO = -i, then the i-th argument had an illegal value |
*> > 0 : DGESVJ did not converge in the maximal allowed number (30) |
*> > 0: DGESVJ did not converge in the maximal allowed number (30) |
*> of sweeps. The output may still be useful. See the |
*> of sweeps. The output may still be useful. See the |
*> description of WORK. |
*> description of WORK. |
*> \endverbatim |
*> \endverbatim |
Line 253
|
Line 255
|
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
|
*> \date September 2012 |
|
* |
* |
*> \ingroup doubleGEcomputational |
*> \ingroup doubleGEcomputational |
* |
* |
Line 279
|
Line 279
|
*> spectral condition number. The best performance of this Jacobi SVD |
*> spectral condition number. The best performance of this Jacobi SVD |
*> procedure is achieved if used in an accelerated version of Drmac and |
*> 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]. |
*> Veselic [5,6], and it is the kernel routine in the SIGMA library [7]. |
*> Some tunning parameters (marked with [TP]) are available for the |
*> Some tuning parameters (marked with [TP]) are available for the |
*> implementer. |
*> implementer. |
*> The computational range for the nonzero singular values is the machine |
*> The computational range for the nonzero singular values is the machine |
*> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even |
*> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even |
Line 335
|
Line 335
|
SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
$ LDV, WORK, LWORK, INFO ) |
$ LDV, WORK, LWORK, INFO ) |
* |
* |
* -- LAPACK computational routine (version 3.4.2) -- |
* -- LAPACK computational routine -- |
* -- 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 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER INFO, LDA, LDV, LWORK, M, MV, N |
INTEGER INFO, LDA, LDV, LWORK, M, MV, N |
Line 374
|
Line 373
|
DOUBLE PRECISION FASTR( 5 ) |
DOUBLE PRECISION FASTR( 5 ) |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC DABS, DMAX1, DMIN1, DBLE, MIN0, DSIGN, DSQRT |
INTRINSIC DABS, MAX, MIN, DBLE, DSIGN, DSQRT |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
* .. |
* .. |
Line 428
|
Line 427
|
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 |
ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN |
ELSE IF( LWORK.LT.MAX( M+N, 6 ) ) THEN |
INFO = -13 |
INFO = -13 |
ELSE |
ELSE |
INFO = 0 |
INFO = 0 |
Line 594
|
Line 593
|
AAPP = ZERO |
AAPP = ZERO |
AAQQ = BIG |
AAQQ = BIG |
DO 4781 p = 1, N |
DO 4781 p = 1, N |
IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) ) |
IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) ) |
AAPP = DMAX1( AAPP, SVA( p ) ) |
AAPP = MAX( AAPP, SVA( p ) ) |
4781 CONTINUE |
4781 CONTINUE |
* |
* |
* #:) Quick return for zero matrix |
* #:) Quick return for zero matrix |
Line 636
|
Line 635
|
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 = MIN( BIG, TEMP1 / AAPP ) |
* AAQQ = AAQQ*TEMP1 |
* AAQQ = AAQQ*TEMP1 |
* AAPP = AAPP*TEMP1 |
* AAPP = AAPP*TEMP1 |
ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN |
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 |
* AAQQ = AAQQ*TEMP1 |
* AAPP = AAPP*TEMP1 |
* AAPP = AAPP*TEMP1 |
ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN |
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 |
* AAQQ = AAQQ*TEMP1 |
* AAPP = AAPP*TEMP1 |
* AAPP = AAPP*TEMP1 |
ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN |
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 |
* AAQQ = AAQQ*TEMP1 |
* AAPP = AAPP*TEMP1 |
* AAPP = AAPP*TEMP1 |
ELSE |
ELSE |
Line 689
|
Line 688
|
* The boundaries are determined dynamically, based on the number of |
* The boundaries are determined dynamically, based on the number of |
* pivots above a threshold. |
* 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 |
*[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 |
* tiling of the p-q loops of pivot pairs. In general, an optimal |
* value of KBL depends on the matrix dimensions and on the |
* value of KBL depends on the matrix dimensions and on the |
Line 701
|
Line 700
|
BLSKIP = KBL**2 |
BLSKIP = KBL**2 |
*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL. |
*[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. |
*[TP] ROWSKIP is a tuning parameter. |
* |
* |
LKAHEAD = 1 |
LKAHEAD = 1 |
Line 712
|
Line 711
|
* invokes cubic convergence. Big part of this cycle is done inside |
* invokes cubic convergence. Big part of this cycle is done inside |
* canonical subspaces of dimensions less than M. |
* 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 |
*[TP] The number of partition levels and the actual partition are |
* tuning parameters. |
* tuning parameters. |
N4 = N / 4 |
N4 = N / 4 |
Line 810
|
Line 809
|
* |
* |
igl = ( ibr-1 )*KBL + 1 |
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 |
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 |
* .. de Rijk's pivoting |
* |
* |
Line 863
|
Line 862
|
* |
* |
PSKIPPED = 0 |
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 ) |
AAQQ = SVA( q ) |
* |
* |
Line 902
|
Line 901
|
END IF |
END IF |
END IF |
END IF |
* |
* |
MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) ) |
MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) ) |
* |
* |
* TO rotate or NOT to rotate, THAT is the question ... |
* TO rotate or NOT to rotate, THAT is the question ... |
* |
* |
Line 935
|
Line 934
|
$ 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( MAX( ZERO, |
$ ONE+T*APOAQ*AAPQ ) ) |
$ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
AAPP = AAPP*DSQRT( MAX( ZERO, |
$ ONE-T*AQOAP*AAPQ ) ) |
$ ONE-T*AQOAP*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( T ) ) |
MXSINJ = MAX( MXSINJ, DABS( T ) ) |
* |
* |
ELSE |
ELSE |
* |
* |
Line 951
|
Line 950
|
CS = DSQRT( ONE / ( ONE+T*T ) ) |
CS = DSQRT( ONE / ( ONE+T*T ) ) |
SN = T*CS |
SN = T*CS |
* |
* |
MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) |
MXSINJ = MAX( MXSINJ, DABS( SN ) ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( MAX( ZERO, |
$ ONE+T*APOAQ*AAPQ ) ) |
$ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
AAPP = AAPP*DSQRT( MAX( ZERO, |
$ ONE-T*AQOAP*AAPQ ) ) |
$ ONE-T*AQOAP*AAPQ ) ) |
* |
* |
APOAQ = WORK( p ) / WORK( q ) |
APOAQ = WORK( p ) / WORK( q ) |
Line 1068
|
Line 1067
|
$ 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( MAX( ZERO, |
$ ONE-AAPQ*AAPQ ) ) |
$ ONE-AAPQ*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
MXSINJ = MAX( MXSINJ, SFMIN ) |
END IF |
END IF |
* END IF ROTOK THEN ... ELSE |
* END IF ROTOK THEN ... ELSE |
* |
* |
Line 1136
|
Line 1135
|
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 + MIN( igl+KBL-1, N ) - p |
END IF |
END IF |
* |
* |
2001 CONTINUE |
2001 CONTINUE |
Line 1156
|
Line 1155
|
* doing the block at ( ibr, jbc ) |
* doing the block at ( ibr, jbc ) |
* |
* |
IJBLSK = 0 |
IJBLSK = 0 |
DO 2100 p = igl, MIN0( igl+KBL-1, N ) |
DO 2100 p = igl, MIN( igl+KBL-1, N ) |
* |
* |
AAPP = SVA( p ) |
AAPP = SVA( p ) |
IF( AAPP.GT.ZERO ) THEN |
IF( AAPP.GT.ZERO ) THEN |
* |
* |
PSKIPPED = 0 |
PSKIPPED = 0 |
* |
* |
DO 2200 q = jgl, MIN0( jgl+KBL-1, N ) |
DO 2200 q = jgl, MIN( jgl+KBL-1, N ) |
* |
* |
AAQQ = SVA( q ) |
AAQQ = SVA( q ) |
IF( AAQQ.GT.ZERO ) THEN |
IF( AAQQ.GT.ZERO ) THEN |
Line 1213
|
Line 1212
|
END IF |
END IF |
END IF |
END IF |
* |
* |
MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) ) |
MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) ) |
* |
* |
* TO rotate or NOT to rotate, THAT is the question ... |
* TO rotate or NOT to rotate, THAT is the question ... |
* |
* |
Line 1241
|
Line 1240
|
$ 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( MAX( ZERO, |
$ ONE+T*APOAQ*AAPQ ) ) |
$ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
AAPP = AAPP*DSQRT( MAX( ZERO, |
$ ONE-T*AQOAP*AAPQ ) ) |
$ ONE-T*AQOAP*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, DABS( T ) ) |
MXSINJ = MAX( MXSINJ, DABS( T ) ) |
ELSE |
ELSE |
* |
* |
* .. choose correct signum for THETA and rotate |
* .. choose correct signum for THETA and rotate |
Line 1256
|
Line 1255
|
$ 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 = MAX( MXSINJ, DABS( SN ) ) |
SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, |
SVA( q ) = AAQQ*DSQRT( MAX( ZERO, |
$ ONE+T*APOAQ*AAPQ ) ) |
$ ONE+T*APOAQ*AAPQ ) ) |
AAPP = AAPP*DSQRT( DMAX1( ZERO, |
AAPP = AAPP*DSQRT( MAX( ZERO, |
$ ONE-T*AQOAP*AAPQ ) ) |
$ ONE-T*AQOAP*AAPQ ) ) |
* |
* |
APOAQ = WORK( p ) / WORK( q ) |
APOAQ = WORK( p ) / WORK( q ) |
Line 1376
|
Line 1375
|
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( MAX( ZERO, |
$ ONE-AAPQ*AAPQ ) ) |
$ ONE-AAPQ*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
MXSINJ = MAX( 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 ) |
Line 1394
|
Line 1393
|
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( MAX( ZERO, |
$ ONE-AAPQ*AAPQ ) ) |
$ ONE-AAPQ*AAPQ ) ) |
MXSINJ = DMAX1( MXSINJ, SFMIN ) |
MXSINJ = MAX( MXSINJ, SFMIN ) |
END IF |
END IF |
END IF |
END IF |
* END IF ROTOK THEN ... ELSE |
* END IF ROTOK THEN ... ELSE |
Line 1466
|
Line 1465
|
ELSE |
ELSE |
* |
* |
IF( AAPP.EQ.ZERO )NOTROT = NOTROT + |
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 |
IF( AAPP.LT.ZERO )NOTROT = 0 |
* |
* |
END IF |
END IF |
Line 1477
|
Line 1476
|
* end of the jbc-loop |
* end of the jbc-loop |
2011 CONTINUE |
2011 CONTINUE |
*2011 bailed out of the jbc-loop |
*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 ) ) |
SVA( p ) = DABS( SVA( p ) ) |
2012 CONTINUE |
2012 CONTINUE |
*** |
*** |