aboutsummaryrefslogtreecommitdiff
path: root/SRC/sgesvj.f
diff options
context:
space:
mode:
authorjason <jason@8a072113-8704-0410-8d35-dd094bca7971>2008-12-30 21:27:12 +0000
committerjason <jason@8a072113-8704-0410-8d35-dd094bca7971>2008-12-30 21:27:12 +0000
commit61e82a389d6bdbdb610d200153937c38c6f6051a (patch)
tree2c16705f27176c833e640589fae73d93af3ffb0b /SRC/sgesvj.f
parentff981f106bde4ce6a74aa4f4a572c943f5a395b2 (diff)
Merged revisions 609-614 via svnmerge from
https://jason@icl.cs.utk.edu/svn/lapack-dev/lapack/branches/SC08-release ........ r609 | julie | 2008-12-16 17:17:52 -0500 (Tue, 16 Dec 2008) | 1 line Polish routines to fit the LAPACK framework and allow manpages generation ........ r610 | langou | 2008-12-19 12:12:38 -0500 (Fri, 19 Dec 2008) | 30 lines bug reported on the forum https://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=854 the complete thread is available at http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/635192e11beadb93# Tobias Burnus also sent us an email: > Hello, > > this was reported at > http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/635192e11beadb93# > > The problem is the line 47: > > 47: IF( M.EQ.0 .OR. A(M, 1).NE.ZERO .OR. A(M, N).NE.ZERO ) THEN > > If M == 0 the one accesses A(0,1) which is invalid as the lower bound is 1 > and not 0. > > Note: Contrary to C there is no left-to-right evaluation of expressions in > Fortran; the order is left to the compiler. One might assume that a smart > compiler does not evaluate "A(M,1)" if "M==0", however, there is nothing in > the standard guarantees this. > > If bounds checks are turned on (see post at the URL above), gfortran aborts > with an out-of-bounds error. ........ r611 | julie | 2008-12-19 15:00:58 -0500 (Fri, 19 Dec 2008) | 5 lines Modify the formatting of the comments. Replace Note and Notes section by Further Details This allow the manpages to be generated corectly. ........ r612 | julie | 2008-12-19 16:29:21 -0500 (Fri, 19 Dec 2008) | 3 lines Reformat the xblas routines comments to be able to generate the manpages ........ r613 | julie | 2008-12-19 16:30:31 -0500 (Fri, 19 Dec 2008) | 1 line Update version number ........ r614 | jason | 2008-12-27 09:44:45 -0500 (Sat, 27 Dec 2008) | 13 lines Fix non-short-circuited tests in ILAxL{C,R}. Fortran doesn't short-circuit logical operators, so the check that the leading dimension /= 0 may not prevent indexing into a 0-length array. Reported by "hes selex" in http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/635192e11beadb93 and forwarded to the LAPACK maintainers by Tobias Burnus <burnus@net-b.de>. Chalk up more bugs found by gfortran's diagnostics! Signed-off-by: Jason Riedy <ejr@cs.berkeley.edu> Cc: Tobias Burnus <burnus@net-b.de> ........
Diffstat (limited to 'SRC/sgesvj.f')
-rw-r--r--SRC/sgesvj.f1576
1 files changed, 867 insertions, 709 deletions
diff --git a/SRC/sgesvj.f b/SRC/sgesvj.f
index 71193ee1..197c4038 100644
--- a/SRC/sgesvj.f
+++ b/SRC/sgesvj.f
@@ -1,5 +1,5 @@
- SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA,
- & MV, V, LDV, WORK, LWORK, INFO )
+ SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
+ + LDV, WORK, LWORK, INFO )
*
* -- LAPACK routine (version 3.2) --
*
@@ -15,19 +15,20 @@
* 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.
*
-* -#- Scalar Arguments -#-
-*
- IMPLICIT NONE
- INTEGER INFO, LDA, LDV, LWORK, M, MV, N
- CHARACTER*1 JOBA, JOBU, JOBV
-*
-* -#- Array Arguments -#-
-*
- REAL A( LDA, * ), SVA( N ), V( LDV, * ), WORK( LWORK )
+ IMPLICIT NONE
+* ..
+* .. Scalar Arguments ..
+ INTEGER INFO, LDA, LDV, LWORK, M, MV, N
+ CHARACTER*1 JOBA, JOBU, JOBV
+* ..
+* .. Array Arguments ..
+ REAL A( LDA, * ), SVA( N ), V( LDV, * ),
+ + WORK( LWORK )
* ..
*
* Purpose
-* ~~~~~~~
+* =======
+*
* SGESVJ 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]
@@ -90,7 +91,7 @@
* drmac@math.hr. Thank you.
*
* Arguments
-* ~~~~~~~~~
+* =========
*
* JOBA (input) CHARACTER* 1
* Specifies the structure of A.
@@ -101,7 +102,6 @@
* 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.
@@ -143,9 +143,7 @@
* On entry, the M-by-N matrix A.
* On exit,
* If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C':
-* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-* If INFO .EQ. 0,
-* ~~~~~~~~~~~~~~~
+* 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
@@ -158,7 +156,6 @@
* TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
* see the description of JOBU.
* If INFO .GT. 0,
-* ~~~~~~~~~~~~~~~
* the procedure SGESVJ did not converge in the given number
* of iterations (sweeps). In that case, the computed
* columns of U may not be orthogonal up to TOL. The output
@@ -166,11 +163,8 @@
* values in SVA(1:N)) and V is still a decomposition of the
* input matrix A in the sense that the residual
* ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
-*
* If JOBU .EQ. 'N':
-* ~~~~~~~~~~~~~~~~~
-* If INFO .EQ. 0
-* ~~~~~~~~~~~~~~
+* If 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
@@ -179,8 +173,7 @@
* numerically orthogonal up to approximately M*EPS. Thus,
* on exit, A contains the columns of U scaled with the
* corresponding singular values.
-* If INFO .GT. 0,
-* ~~~~~~~~~~~~~~~
+* If INFO .GT. 0 :
* the procedure SGESVJ did not converge in the given number
* of iterations (sweeps).
*
@@ -189,22 +182,18 @@
*
* SVA (workspace/output) REAL array, dimension (N)
* On exit,
-* If INFO .EQ. 0,
-* ~~~~~~~~~~~~~~~
+* 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,
-* ~~~~~~~~~~~~~~~
+* If INFO .GT. 0 :
* the procedure SGESVJ did not converge in the given number of
* iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
*
@@ -227,8 +216,7 @@
*
* WORK (input/workspace/output) REAL array, dimension max(4,M+N).
* On entry,
-* If JOBU .EQ. 'C',
-* ~~~~~~~~~~~~~~~~~
+* 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=SLAMCH('E').
@@ -261,55 +249,55 @@
* > 0 : SGESVJ did not converge in the maximal allowed number (30)
* of sweeps. The output may still be useful. See the
* description of WORK.
+* =====================================================================
+*
+* .. Local Parameters ..
+ REAL ZERO, HALF, ONE, TWO
+ PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
+ + TWO = 2.0E0 )
+ INTEGER NSWEEP
+ PARAMETER ( NSWEEP = 30 )
+* ..
+* .. Local Scalars ..
+ REAL 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
+ 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
+ LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
+ + RSVEC, UCTOL, UPPER
+* ..
+* .. Local Arrays ..
+ REAL FASTR( 5 )
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, AMAX1, AMIN1, FLOAT, MIN0, SIGN, SQRT
+* ..
+* .. External Functions ..
+* from BLAS
+ REAL SDOT, SNRM2
+ EXTERNAL SDOT, SNRM2
+ INTEGER ISAMAX
+ EXTERNAL ISAMAX
+* from LAPACK
+ REAL SLAMCH
+ EXTERNAL SLAMCH
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+* from BLAS
+ EXTERNAL SAXPY, SCOPY, SROTM, SSCAL, SSWAP
+* from LAPACK
+ EXTERNAL SLASCL, SLASET, SLASSQ, XERBLA
*
-* Local Parameters
-*
- REAL ZERO, HALF, ONE, TWO
- PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0, TWO = 2.0E0 )
- INTEGER NSWEEP
- PARAMETER ( NSWEEP = 30 )
-*
-* Local Scalars
-*
- REAL 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
- 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
- LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
- & RSVEC, UCTOL, UPPER
-*
-* Local Arrays
-*
- REAL FASTR(5)
-*
-* Intrinsic Functions
-*
- INTRINSIC ABS, AMAX1, AMIN1, FLOAT, MIN0, SIGN, SQRT
-*
-* External Functions
-* .. from BLAS
- REAL SDOT, SNRM2
- EXTERNAL SDOT, SNRM2
- INTEGER ISAMAX
- EXTERNAL ISAMAX
-* .. from LAPACK
- REAL SLAMCH
- EXTERNAL SLAMCH
- LOGICAL LSAME
- EXTERNAL LSAME
-*
-* External Subroutines
-* .. from BLAS
- EXTERNAL SAXPY, SCOPY, SROTM, SSCAL, SSWAP
-* .. from LAPACK
- EXTERNAL SLASCL, SLASET, SLASSQ, XERBLA
-*
- EXTERNAL SGSVJ0, SGSVJ1
+ EXTERNAL SGSVJ0, SGSVJ1
+* ..
+* .. Executable Statements ..
*
* Test the input arguments
*
@@ -320,40 +308,40 @@
UPPER = LSAME( JOBA, 'U' )
LOWER = LSAME( JOBA, 'L' )
*
- IF ( .NOT.( UPPER .OR. LOWER .OR. LSAME(JOBA,'G') ) ) THEN
- INFO = - 1
- ELSE IF ( .NOT.( LSVEC .OR. UCTOL .OR. LSAME(JOBU,'N') ) ) THEN
- INFO = - 2
- ELSE IF ( .NOT.( RSVEC .OR. APPLV .OR. LSAME(JOBV,'N') ) ) THEN
- INFO = - 3
- ELSE IF ( M .LT. 0 ) THEN
- INFO = - 4
- ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
- INFO = - 5
- ELSE IF ( LDA .LT. M ) THEN
- INFO = - 7
- ELSE IF ( MV .LT. 0 ) THEN
- INFO = - 9
- ELSE IF ( ( RSVEC .AND. (LDV .LT. N ) ) .OR.
- & ( APPLV .AND. (LDV .LT. MV) ) ) THEN
+ IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
+ INFO = -2
+ ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
+ INFO = -3
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
+ INFO = -5
+ ELSE IF( LDA.LT.M ) THEN
+ INFO = -7
+ ELSE IF( MV.LT.0 ) THEN
+ INFO = -9
+ ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
+ + ( 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
- INFO = - 13
+ ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
+ INFO = -12
+ ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN
+ INFO = -13
ELSE
- INFO = 0
+ INFO = 0
END IF
*
* #:(
- IF ( INFO .NE. 0 ) THEN
+ IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SGESVJ', -INFO )
RETURN
END IF
*
* #:) Quick return for void matrix
*
- IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
+ IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN
*
* Set numerical parameters
* The stopping criterion for Jacobi rotations is
@@ -362,45 +350,45 @@
*
* where EPS is the round-off and CTOL is defined as follows:
*
- IF ( UCTOL ) THEN
+ IF( UCTOL ) THEN
* ... user controlled
- CTOL = WORK(1)
+ CTOL = WORK( 1 )
ELSE
* ... default
- IF ( LSVEC .OR. RSVEC .OR. APPLV ) THEN
- CTOL = SQRT(FLOAT(M))
+ IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
+ CTOL = SQRT( FLOAT( M ) )
ELSE
- CTOL = FLOAT(M)
+ CTOL = FLOAT( M )
END IF
END IF
* ... and the machine dependent parameters are
*[!] (Make sure that SLAMCH() works properly on the target machine.)
*
- EPSILON = SLAMCH('Epsilon')
- ROOTEPS = SQRT(EPSILON)
- SFMIN = SLAMCH('SafeMinimum')
- ROOTSFMIN = SQRT(SFMIN)
- SMALL = SFMIN / EPSILON
- BIG = SLAMCH('Overflow')
- ROOTBIG = ONE / ROOTSFMIN
- LARGE = BIG / SQRT(FLOAT(M*N))
- BIGTHETA = ONE / ROOTEPS
-*
- TOL = CTOL * EPSILON
- ROOTTOL = SQRT(TOL)
-*
- IF ( FLOAT(M)*EPSILON .GE. ONE ) THEN
- INFO = - 5
+ EPSILON = SLAMCH( 'Epsilon' )
+ ROOTEPS = SQRT( EPSILON )
+ SFMIN = SLAMCH( 'SafeMinimum' )
+ ROOTSFMIN = SQRT( SFMIN )
+ SMALL = SFMIN / EPSILON
+ BIG = SLAMCH( 'Overflow' )
+ ROOTBIG = ONE / ROOTSFMIN
+ LARGE = BIG / SQRT( FLOAT( M*N ) )
+ BIGTHETA = ONE / ROOTEPS
+*
+ TOL = CTOL*EPSILON
+ ROOTTOL = SQRT( TOL )
+*
+ IF( FLOAT( M )*EPSILON.GE.ONE ) THEN
+ INFO = -5
CALL XERBLA( 'SGESVJ', -INFO )
RETURN
END IF
*
* Initialize the right singular vector matrix.
*
- IF ( RSVEC ) THEN
+ IF( RSVEC ) THEN
MVL = N
CALL SLASET( 'A', MVL, N, ZERO, ONE, V, LDV )
- ELSE IF ( APPLV ) THEN
+ ELSE IF( APPLV ) THEN
MVL = MV
END IF
RSVEC = RSVEC .OR. APPLV
@@ -414,56 +402,56 @@
* SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
* in A are detected, the procedure returns with INFO=-6.
*
- SCALE = ONE / SQRT(FLOAT(M)*FLOAT(N))
- NOSCALE = .TRUE.
- GOSCALE = .TRUE.
+ SCALE = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
+ NOSCALE = .TRUE.
+ GOSCALE = .TRUE.
*
- IF ( LOWER ) THEN
+ IF( LOWER ) THEN
* the input matrix is M-by-N lower triangular (trapezoidal)
DO 1874 p = 1, N
AAPP = ZERO
AAQQ = ZERO
- CALL SLASSQ( M-p+1, A(p,p), 1, AAPP, AAQQ )
- IF ( AAPP .GT. BIG ) THEN
- INFO = - 6
+ CALL SLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
+ IF( AAPP.GT.BIG ) THEN
+ INFO = -6
CALL XERBLA( 'SGESVJ', -INFO )
RETURN
END IF
- AAQQ = SQRT(AAQQ)
- IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCALE ) THEN
- SVA(p) = AAPP * AAQQ
+ AAQQ = SQRT( AAQQ )
+ IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
+ SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA(p) = AAPP * ( AAQQ * SCALE )
- IF ( GOSCALE ) THEN
+ SVA( p ) = AAPP*( AAQQ*SCALE )
+ IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 1873 q = 1, p - 1
- SVA(q) = SVA(q)*SCALE
+ SVA( q ) = SVA( q )*SCALE
1873 CONTINUE
END IF
END IF
1874 CONTINUE
- ELSE IF ( UPPER ) THEN
+ ELSE IF( UPPER ) THEN
* the input matrix is M-by-N upper triangular (trapezoidal)
DO 2874 p = 1, N
AAPP = ZERO
AAQQ = ZERO
- CALL SLASSQ( p, A(1,p), 1, AAPP, AAQQ )
- IF ( AAPP .GT. BIG ) THEN
- INFO = - 6
+ CALL SLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
+ IF( AAPP.GT.BIG ) THEN
+ INFO = -6
CALL XERBLA( 'SGESVJ', -INFO )
RETURN
END IF
- AAQQ = SQRT(AAQQ)
- IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCALE ) THEN
- SVA(p) = AAPP * AAQQ
+ AAQQ = SQRT( AAQQ )
+ IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
+ SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA(p) = AAPP * ( AAQQ * SCALE )
- IF ( GOSCALE ) THEN
+ SVA( p ) = AAPP*( AAQQ*SCALE )
+ IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 2873 q = 1, p - 1
- SVA(q) = SVA(q)*SCALE
+ SVA( q ) = SVA( q )*SCALE
2873 CONTINUE
END IF
END IF
@@ -473,29 +461,29 @@
DO 3874 p = 1, N
AAPP = ZERO
AAQQ = ZERO
- CALL SLASSQ( M, A(1,p), 1, AAPP, AAQQ )
- IF ( AAPP .GT. BIG ) THEN
- INFO = - 6
+ CALL SLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
+ IF( AAPP.GT.BIG ) THEN
+ INFO = -6
CALL XERBLA( 'SGESVJ', -INFO )
RETURN
END IF
- AAQQ = SQRT(AAQQ)
- IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCALE ) THEN
- SVA(p) = AAPP * AAQQ
+ AAQQ = SQRT( AAQQ )
+ IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
+ SVA( p ) = AAPP*AAQQ
ELSE
NOSCALE = .FALSE.
- SVA(p) = AAPP * ( AAQQ * SCALE )
- IF ( GOSCALE ) THEN
+ SVA( p ) = AAPP*( AAQQ*SCALE )
+ IF( GOSCALE ) THEN
GOSCALE = .FALSE.
DO 3873 q = 1, p - 1
- SVA(q) = SVA(q)*SCALE
+ SVA( q ) = SVA( q )*SCALE
3873 CONTINUE
END IF
END IF
3874 CONTINUE
END IF
*
- IF ( NOSCALE ) SCALE = ONE
+ IF( NOSCALE )SCALE = ONE
*
* Move the smaller part of the spectrum from the underflow threshold
*(!) Start by determining the position of the nonzero entries of the
@@ -504,61 +492,61 @@
AAPP = ZERO
AAQQ = BIG
DO 4781 p = 1, N
- IF ( SVA(p) .NE. ZERO ) AAQQ = AMIN1( AAQQ, SVA(p) )
- AAPP = AMAX1( AAPP, SVA(p) )
+ IF( SVA( p ).NE.ZERO )AAQQ = AMIN1( AAQQ, SVA( p ) )
+ AAPP = AMAX1( AAPP, SVA( p ) )
4781 CONTINUE
*
* #:) Quick return for zero matrix
*
- IF ( AAPP .EQ. ZERO ) THEN
- IF ( LSVEC ) CALL SLASET( 'G', M, N, ZERO, ONE, A, LDA )
- WORK(1) = ONE
- WORK(2) = ZERO
- WORK(3) = ZERO
- WORK(4) = ZERO
- WORK(5) = ZERO
- WORK(6) = ZERO
+ IF( AAPP.EQ.ZERO ) THEN
+ IF( LSVEC )CALL SLASET( 'G', M, N, ZERO, ONE, A, LDA )
+ WORK( 1 ) = ONE
+ WORK( 2 ) = ZERO
+ WORK( 3 ) = ZERO
+ WORK( 4 ) = ZERO
+ WORK( 5 ) = ZERO
+ WORK( 6 ) = ZERO
RETURN
END IF
*
* #:) Quick return for one-column matrix
*
- IF ( N .EQ. 1 ) THEN
- IF ( LSVEC )
- & CALL SLASCL( 'G',0,0,SVA(1),SCALE,M,1,A(1,1),LDA,IERR )
- WORK(1) = ONE / SCALE
- IF ( SVA(1) .GE. SFMIN ) THEN
- WORK(2) = ONE
+ IF( N.EQ.1 ) THEN
+ IF( LSVEC )CALL SLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,
+ + A( 1, 1 ), LDA, IERR )
+ WORK( 1 ) = ONE / SCALE
+ IF( SVA( 1 ).GE.SFMIN ) THEN
+ WORK( 2 ) = ONE
ELSE
- WORK(2) = ZERO
+ WORK( 2 ) = ZERO
END IF
- WORK(3) = ZERO
- WORK(4) = ZERO
- WORK(5) = ZERO
- WORK(6) = ZERO
+ WORK( 3 ) = ZERO
+ WORK( 4 ) = ZERO
+ WORK( 5 ) = ZERO
+ WORK( 6 ) = ZERO
RETURN
END IF
*
* Protect small singular values from underflow, and try to
* avoid underflows/overflows in computing Jacobi rotations.
*
- SN = SQRT( SFMIN / EPSILON )
- TEMP1 = SQRT( BIG / FLOAT(N) )
- IF ( (AAPP.LE.SN).OR.(AAQQ.GE.TEMP1)
- & .OR.((SN.LE.AAQQ).AND.(AAPP.LE.TEMP1)) ) THEN
- TEMP1 = AMIN1(BIG,TEMP1/AAPP)
+ SN = SQRT( SFMIN / EPSILON )
+ TEMP1 = SQRT( BIG / FLOAT( N ) )
+ IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
+ + ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
+ TEMP1 = AMIN1( BIG, TEMP1 / AAPP )
* AAQQ = AAQQ*TEMP1
* AAPP = AAPP*TEMP1
- ELSE IF ( (AAQQ.LE.SN).AND.(AAPP.LE.TEMP1) ) THEN
- TEMP1 = AMIN1( SN / AAQQ, BIG/(AAPP*SQRT(FLOAT(N))) )
+ ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
+ TEMP1 = AMIN1( SN / AAQQ, BIG / ( AAPP*SQRT( FLOAT( N ) ) ) )
* AAQQ = AAQQ*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 = AMAX1( SN / AAQQ, TEMP1 / AAPP )
* AAQQ = AAQQ*TEMP1
* AAPP = AAPP*TEMP1
- ELSE IF ( (AAQQ.LE.SN).AND.(AAPP.GE.TEMP1) ) THEN
- TEMP1 = AMIN1( SN / AAQQ, BIG / (SQRT(FLOAT(N))*AAPP))
+ ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
+ TEMP1 = AMIN1( SN / AAQQ, BIG / ( SQRT( FLOAT( N ) )*AAPP ) )
* AAQQ = AAQQ*TEMP1
* AAPP = AAPP*TEMP1
ELSE
@@ -567,27 +555,27 @@
*
* Scale, if necessary
*
- IF ( TEMP1 .NE. ONE ) THEN
+ IF( TEMP1.NE.ONE ) THEN
CALL SLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
END IF
- SCALE = TEMP1 * SCALE
- IF ( SCALE .NE. ONE ) THEN
+ SCALE = TEMP1*SCALE
+ IF( SCALE.NE.ONE ) THEN
CALL SLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )
SCALE = ONE / SCALE
END IF
*
* Row-cyclic Jacobi SVD algorithm with column pivoting
*
- EMPTSW = ( N * ( N - 1 ) ) / 2
- NOTROT = 0
- FASTR(1) = ZERO
+ EMPTSW = ( N*( N-1 ) ) / 2
+ NOTROT = 0
+ FASTR( 1 ) = ZERO
*
* A is represented in factored form A = A * diag(WORK), where diag(WORK)
* is initialized to identity. WORK is updated during fast scaled
* rotations.
*
DO 1868 q = 1, N
- WORK(q) = ONE
+ WORK( q ) = ONE
1868 CONTINUE
*
*
@@ -606,7 +594,7 @@
* parameters of the computer's memory.
*
NBL = N / KBL
- IF ( ( NBL * KBL ) .NE. N ) NBL = NBL + 1
+ IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
*
BLSKIP = KBL**2
*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
@@ -622,19 +610,19 @@
* 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.MAX0( 64, 4*KBL ) ) ) THEN
*[TP] The number of partition levels and the actual partition are
* tuning parameters.
- N4 = N / 4
- N2 = N / 2
- N34 = 3 * N4
- IF ( APPLV ) THEN
- q = 0
- ELSE
- q = 1
- END IF
+ N4 = N / 4
+ N2 = N / 2
+ N34 = 3*N4
+ IF( APPLV ) THEN
+ q = 0
+ ELSE
+ q = 1
+ END IF
*
- IF ( LOWER ) THEN
+ IF( LOWER ) THEN
*
* This works very well on lower triangular matrices, in particular
* in the framework of the preconditioned Jacobi SVD (xGEJSV).
@@ -644,92 +632,103 @@
* [+ + x 0] actually work on [x 0] [x 0]
* [+ + x x] [x x]. [x x]
*
- CALL SGSVJ0(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 )
+ CALL SGSVJ0( 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 )
*
- CALL SGSVJ0( 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 )
+ CALL SGSVJ0( 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 )
*
- CALL SGSVJ1( 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 )
+ CALL SGSVJ1( 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 )
*
- CALL SGSVJ0( 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 )
+ CALL SGSVJ0( 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 )
*
- CALL SGSVJ0( JOBV,M,N4,A,LDA,WORK,SVA,MVL,V,LDV,EPSILON,
- & SFMIN,TOL,1,WORK(N+1),LWORK-N,IERR )
+ CALL SGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
+ + EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
+ + IERR )
*
- CALL SGSVJ1( JOBV,M,N2,N4,A,LDA,WORK,SVA,MVL,V,LDV,EPSILON,
- & SFMIN,TOL,1,WORK(N+1),LWORK-N,IERR )
+ CALL SGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
+ + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+ + LWORK-N, IERR )
*
*
- ELSE IF ( UPPER ) THEN
+ ELSE IF( UPPER ) THEN
*
*
- CALL SGSVJ0( JOBV,N4,N4,A,LDA,WORK,SVA,MVL,V,LDV,EPSILON,
- & SFMIN,TOL,2,WORK(N+1),LWORK-N,IERR )
+ CALL SGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
+ + EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
+ + IERR )
*
- CALL SGSVJ0(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 )
+ CALL SGSVJ0( 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 )
*
- CALL SGSVJ1( JOBV,N2,N2,N4,A,LDA,WORK,SVA,MVL,V,LDV,EPSILON,
- & SFMIN,TOL,1,WORK(N+1),LWORK-N,IERR )
+ CALL SGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
+ + LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
+ + LWORK-N, IERR )
*
- CALL SGSVJ0( 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 )
+ CALL SGSVJ0( 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 )
- END IF
+ END IF
*
END IF
*
-* -#- Row-cyclic pivot strategy with de Rijk's pivoting -#-
+* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
*
DO 1993 i = 1, NSWEEP
* .. go go go ...
*
- MXAAPQ = ZERO
- MXSINJ = ZERO
- ISWROT = 0
+ MXAAPQ = ZERO
+ MXSINJ = ZERO
+ ISWROT = 0
*
- NOTROT = 0
- PSKIPPED = 0
+ NOTROT = 0
+ PSKIPPED = 0
*
* Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
* 1 <= p < q <= N. This is the first step toward a blocked implementation
* of the rotations. New implementation, based on block transformations,
* is under development.
*
- DO 2000 ibr = 1, NBL
+ DO 2000 ibr = 1, NBL
*
- igl = ( ibr - 1 ) * KBL + 1
+ igl = ( ibr-1 )*KBL + 1
*
- DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL - ibr )
+ DO 1002 ir1 = 0, MIN0( 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, MIN0( igl+KBL-1, N-1 )
*
* .. de Rijk's pivoting
*
- q = ISAMAX( N-p+1, SVA(p), 1 ) + p - 1
- IF ( p .NE. q ) THEN
- CALL SSWAP( M, A(1,p), 1, A(1,q), 1 )
- IF ( RSVEC ) CALL SSWAP( MVL, V(1,p), 1, V(1,q), 1 )
- TEMP1 = SVA(p)
- SVA(p) = SVA(q)
- SVA(q) = TEMP1
- TEMP1 = WORK(p)
- WORK(p) = WORK(q)
- WORK(q) = TEMP1
- END IF
-*
- IF ( ir1 .EQ. 0 ) THEN
+ q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
+ IF( p.NE.q ) THEN
+ CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
+ IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1,
+ + V( 1, q ), 1 )
+ TEMP1 = SVA( p )
+ SVA( p ) = SVA( q )
+ SVA( q ) = TEMP1
+ TEMP1 = WORK( p )
+ WORK( p ) = WORK( q )
+ WORK( q ) = TEMP1
+ END IF
+*
+ IF( ir1.EQ.0 ) THEN
*
* Column norms are periodically updated by explicit
* norm computation.
@@ -743,506 +742,665 @@
* If properly implemented SNRM2 is available, the IF-THEN-ELSE
* below should read "AAPP = SNRM2( M, A(1,p), 1 ) * WORK(p)".
*
- IF ((SVA(p) .LT. ROOTBIG) .AND. (SVA(p) .GT. ROOTSFMIN)) THEN
- SVA(p) = SNRM2( M, A(1,p), 1 ) * WORK(p)
- ELSE
- TEMP1 = ZERO
- AAPP = ZERO
- CALL SLASSQ( M, A(1,p), 1, TEMP1, AAPP )
- SVA(p) = TEMP1 * SQRT(AAPP) * WORK(p)
- END IF
- AAPP = SVA(p)
- ELSE
- AAPP = SVA(p)
- END IF
-*
- IF ( AAPP .GT. ZERO ) THEN
-*
- PSKIPPED = 0
-*
- DO 2002 q = p + 1, MIN0( igl + KBL - 1, N )
-*
- AAQQ = SVA(q)
-*
- IF ( AAQQ .GT. ZERO ) THEN
-*
- AAPP0 = AAPP
- IF ( AAQQ .GE. ONE ) THEN
- ROTOK = ( SMALL*AAPP ) .LE. AAQQ
- IF ( AAPP .LT. ( BIG / AAQQ ) ) THEN
- AAPQ = ( SDOT(M, A(1,p), 1, A(1,q), 1 ) *
- & WORK(p) * WORK(q) / AAQQ ) / AAPP
- ELSE
- CALL SCOPY( M, A(1,p), 1, WORK(N+1), 1 )
- CALL SLASCL( 'G', 0, 0, AAPP, WORK(p), M,
- & 1, WORK(N+1), LDA, IERR )
- AAPQ = SDOT( M, WORK(N+1),1, A(1,q),1 )*WORK(q) / AAQQ
- END IF
- ELSE
- ROTOK = AAPP .LE. ( AAQQ / SMALL )
- IF ( AAPP .GT. ( SMALL / AAQQ ) ) THEN
- AAPQ = ( SDOT( M, A(1,p), 1, A(1,q), 1 ) *
- & WORK(p) * WORK(q) / AAQQ ) / AAPP
- ELSE
- CALL SCOPY( M, A(1,q), 1, WORK(N+1), 1 )
- CALL SLASCL( 'G', 0, 0, AAQQ, WORK(q), M,
- & 1, WORK(N+1), LDA, IERR )
- AAPQ = SDOT( M, WORK(N+1),1, A(1,p),1 )*WORK(p) / AAPP
- END IF
- END IF
+ IF( ( SVA( p ).LT.ROOTBIG ) .AND.
+ + ( SVA( p ).GT.ROOTSFMIN ) ) THEN
+ SVA( p ) = SNRM2( M, A( 1, p ), 1 )*WORK( p )
+ ELSE
+ TEMP1 = ZERO
+ AAPP = ZERO
+ CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
+ SVA( p ) = TEMP1*SQRT( AAPP )*WORK( p )
+ END IF
+ AAPP = SVA( p )
+ ELSE
+ AAPP = SVA( p )
+ END IF
+*
+ IF( AAPP.GT.ZERO ) THEN
+*
+ PSKIPPED = 0
+*
+ DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
+*
+ AAQQ = SVA( q )
+*
+ IF( AAQQ.GT.ZERO ) THEN
+*
+ AAPP0 = AAPP
+ IF( AAQQ.GE.ONE ) THEN
+ ROTOK = ( SMALL*AAPP ).LE.AAQQ
+ IF( AAPP.LT.( BIG / AAQQ ) ) THEN
+ AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
+ + q ), 1 )*WORK( p )*WORK( q ) /
+ + AAQQ ) / AAPP
+ ELSE
+ CALL SCOPY( M, A( 1, p ), 1,
+ + WORK( N+1 ), 1 )
+ CALL SLASCL( 'G', 0, 0, AAPP,
+ + WORK( p ), M, 1,
+ + WORK( N+1 ), LDA, IERR )
+ AAPQ = SDOT( M, WORK( N+1 ), 1,
+ + A( 1, q ), 1 )*WORK( q ) / AAQQ
+ END IF
+ ELSE
+ ROTOK = AAPP.LE.( AAQQ / SMALL )
+ IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
+ AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
+ + q ), 1 )*WORK( p )*WORK( q ) /
+ + AAQQ ) / AAPP
+ ELSE
+ CALL SCOPY( M, A( 1, q ), 1,
+ + WORK( N+1 ), 1 )
+ CALL SLASCL( 'G', 0, 0, AAQQ,
+ + WORK( q ), M, 1,
+ + WORK( N+1 ), LDA, IERR )
+ AAPQ = SDOT( M, WORK( N+1 ), 1,
+ + A( 1, p ), 1 )*WORK( p ) / AAPP
+ END IF
+ END IF
*
- MXAAPQ = AMAX1( MXAAPQ, ABS(AAPQ) )
+ MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
*
* TO rotate or NOT to rotate, THAT is the question ...
*
- IF ( ABS( AAPQ ) .GT. TOL ) THEN
+ IF( ABS( AAPQ ).GT.TOL ) THEN
*
* .. rotate
*[RTD] ROTATED = ROTATED + ONE
*
- IF ( ir1 .EQ. 0 ) THEN
- NOTROT = 0
- PSKIPPED = 0
- ISWROT = ISWROT + 1
- END IF
-*
- IF ( ROTOK ) THEN
-*
- AQOAP = AAQQ / AAPP
- APOAQ = AAPP / AAQQ
- THETA = - HALF * ABS( AQOAP - APOAQ ) / AAPQ
-*
- IF ( ABS( THETA ) .GT. BIGTHETA ) THEN
-*
- T = HALF / THETA
- FASTR(3) = T * WORK(p) / WORK(q)
- FASTR(4) = - T * WORK(q) / WORK(p)
- CALL SROTM( M, A(1,p), 1, A(1,q), 1, FASTR )
- IF ( RSVEC )
- & CALL SROTM( MVL, V(1,p), 1, V(1,q), 1, FASTR )
- SVA(q) = AAQQ*SQRT( AMAX1(ZERO,ONE + T*APOAQ*AAPQ) )
- AAPP = AAPP*SQRT( ONE - T*AQOAP*AAPQ )
- MXSINJ = AMAX1( MXSINJ, ABS(T) )
-*
- ELSE
+ IF( ir1.EQ.0 ) THEN
+ NOTROT = 0
+ PSKIPPED = 0
+ ISWROT = ISWROT + 1
+ END IF
+*
+ IF( ROTOK ) THEN
+*
+ AQOAP = AAQQ / AAPP
+ APOAQ = AAPP / AAQQ
+ THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
+*
+ IF( ABS( THETA ).GT.BIGTHETA ) THEN
+*
+ T = HALF / THETA
+ FASTR( 3 ) = T*WORK( p ) / WORK( q )
+ FASTR( 4 ) = -T*WORK( q ) /
+ + WORK( p )
+ CALL SROTM( M, A( 1, p ), 1,
+ + A( 1, q ), 1, FASTR )
+ IF( RSVEC )CALL SROTM( MVL,
+ + V( 1, p ), 1,
+ + V( 1, q ), 1,
+ + FASTR )
+ SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ + ONE+T*APOAQ*AAPQ ) )
+ AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+ MXSINJ = AMAX1( MXSINJ, ABS( T ) )
+*
+ ELSE
*
* .. choose correct signum for THETA and rotate
*
- THSIGN = - SIGN(ONE,AAPQ)
- T = ONE / ( THETA + THSIGN*SQRT(ONE+THETA*THETA) )
- CS = SQRT( ONE / ( ONE + T*T ) )
- SN = T * CS
-*
- MXSINJ = AMAX1( MXSINJ, ABS(SN) )
- SVA(q) = AAQQ*SQRT( AMAX1(ZERO, ONE+T*APOAQ*AAPQ) )
- AAPP = AAPP*SQRT( AMAX1(ZERO, ONE-T*AQOAP*AAPQ) )
-*
- APOAQ = WORK(p) / WORK(q)
- AQOAP = WORK(q) / WORK(p)
- IF ( WORK(p) .GE. ONE ) THEN
- IF ( WORK(q) .GE. ONE ) THEN
- FASTR(3) = T * APOAQ
- FASTR(4) = - T * AQOAP
- WORK(p) = WORK(p) * CS
- WORK(q) = WORK(q) * CS
- CALL SROTM( M, A(1,p),1, A(1,q),1, FASTR )
- IF ( RSVEC )
- & CALL SROTM( MVL, V(1,p),1, V(1,q),1, FASTR )
- ELSE
- CALL SAXPY( M, -T*AQOAP, A(1,q),1, A(1,p),1 )
- CALL SAXPY( M, CS*SN*APOAQ, A(1,p),1, A(1,q),1 )
- WORK(p) = WORK(p) * CS
- WORK(q) = WORK(q) / CS
- IF ( RSVEC ) THEN
- CALL SAXPY(MVL, -T*AQOAP, V(1,q),1,V(1,p),1)
- CALL SAXPY(MVL,CS*SN*APOAQ, V(1,p),1,V(1,q),1)
- END IF
- END IF
- ELSE
- IF ( WORK(q) .GE. ONE ) THEN
- CALL SAXPY( M, T*APOAQ, A(1,p),1, A(1,q),1 )
- CALL SAXPY( M,-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 SAXPY(MVL, T*APOAQ, V(1,p),1,V(1,q),1)
- CALL SAXPY(MVL,-CS*SN*AQOAP,V(1,q),1,V(1,p),1)
- END IF
- ELSE
- IF ( WORK(p) .GE. WORK(q) ) THEN
- CALL SAXPY( M,-T*AQOAP, A(1,q),1,A(1,p),1 )
- CALL SAXPY( M,CS*SN*APOAQ,A(1,p),1,A(1,q),1 )
- WORK(p) = WORK(p) * CS
- WORK(q) = WORK(q) / CS
- IF ( RSVEC ) THEN
- CALL SAXPY(MVL, -T*AQOAP, V(1,q),1,V(1,p),1)
- CALL SAXPY(MVL,CS*SN*APOAQ,V(1,p),1,V(1,q),1)
- END IF
- ELSE
- CALL SAXPY( M, T*APOAQ, A(1,p),1,A(1,q),1)
- CALL SAXPY( M,-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 SAXPY(MVL, T*APOAQ, V(1,p),1,V(1,q),1)
- CALL SAXPY(MVL,-CS*SN*AQOAP,V(1,q),1,V(1,p),1)
- END IF
- END IF
- END IF
- ENDIF
- END IF
-*
- ELSE
+ THSIGN = -SIGN( ONE, AAPQ )
+ T = ONE / ( THETA+THSIGN*
+ + SQRT( ONE+THETA*THETA ) )
+ CS = SQRT( ONE / ( ONE+T*T ) )
+ SN = T*CS
+*
+ MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
+ SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ + ONE+T*APOAQ*AAPQ ) )
+ AAPP = AAPP*SQRT( AMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
+*
+ APOAQ = WORK( p ) / WORK( q )
+ AQOAP = WORK( q ) / WORK( p )
+ IF( WORK( p ).GE.ONE ) THEN
+ IF( WORK( q ).GE.ONE ) THEN
+ FASTR( 3 ) = T*APOAQ
+ FASTR( 4 ) = -T*AQOAP
+ WORK( p ) = WORK( p )*CS
+ WORK( q ) = WORK( q )*CS
+ CALL SROTM( M, A( 1, p ), 1,
+ + A( 1, q ), 1,
+ + FASTR )
+ IF( RSVEC )CALL SROTM( MVL,
+ + V( 1, p ), 1, V( 1, q ),
+ + 1, FASTR )
+ ELSE
+ CALL SAXPY( M, -T*AQOAP,
+ + A( 1, q ), 1,
+ + A( 1, p ), 1 )
+ CALL SAXPY( M, CS*SN*APOAQ,
+ + A( 1, p ), 1,
+ + A( 1, q ), 1 )
+ WORK( p ) = WORK( p )*CS
+ WORK( q ) = WORK( q ) / CS
+ IF( RSVEC ) THEN
+ CALL SAXPY( MVL, -T*AQOAP,
+ + V( 1, q ), 1,
+ + V( 1, p ), 1 )
+ CALL SAXPY( MVL,
+ + CS*SN*APOAQ,
+ + V( 1, p ), 1,
+ + V( 1, q ), 1 )
+ END IF
+ END IF
+ ELSE
+ IF( WORK( q ).GE.ONE ) THEN
+ CALL SAXPY( M, T*APOAQ,
+ + A( 1, p ), 1,
+ + A( 1, q ), 1 )
+ CALL SAXPY( M, -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 SAXPY( MVL, T*APOAQ,
+ + V( 1, p ), 1,
+ + V( 1, q ), 1 )
+ CALL SAXPY( MVL,
+ + -CS*SN*AQOAP,
+ + V( 1, q ), 1,
+ + V( 1, p ), 1 )
+ END IF
+ ELSE
+ IF( WORK( p ).GE.WORK( q ) )
+ + THEN
+ CALL SAXPY( M, -T*AQOAP,
+ + A( 1, q ), 1,
+ + A( 1, p ), 1 )
+ CALL SAXPY( M, CS*SN*APOAQ,
+ + A( 1, p ), 1,
+ + A( 1, q ), 1 )
+ WORK( p ) = WORK( p )*CS
+ WORK( q ) = WORK( q ) / CS
+ IF( RSVEC ) THEN
+ CALL SAXPY( MVL,
+ + -T*AQOAP,
+ + V( 1, q ), 1,
+ + V( 1, p ), 1 )
+ CALL SAXPY( MVL,
+ + CS*SN*APOAQ,
+ + V( 1, p ), 1,
+ + V( 1, q ), 1 )
+ END IF
+ ELSE
+ CALL SAXPY( M, T*APOAQ,
+ + A( 1, p ), 1,
+ + A( 1, q ), 1 )
+ CALL SAXPY( M,
+ + -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 SAXPY( MVL,
+ + T*APOAQ, V( 1, p ),
+ + 1, V( 1, q ), 1 )
+ CALL SAXPY( MVL,
+ + -CS*SN*AQOAP,
+ + V( 1, q ), 1,
+ + V( 1, p ), 1 )
+ END IF
+ END IF
+ END IF
+ END IF
+ END IF
+*
+ ELSE
* .. have to use modified Gram-Schmidt like transformation
- CALL SCOPY( M, A(1,p), 1, WORK(N+1), 1 )
- CALL SLASCL( 'G',0,0,AAPP,ONE,M,1,WORK(N+1),LDA,IERR )
- CALL SLASCL( 'G',0,0,AAQQ,ONE,M,1, A(1,q),LDA,IERR )
- TEMP1 = -AAPQ * WORK(p) / WORK(q)
- CALL SAXPY ( M, TEMP1, WORK(N+1), 1, A(1,q), 1 )
- CALL SLASCL( 'G',0,0,ONE,AAQQ,M,1, A(1,q),LDA,IERR )
- SVA(q) = AAQQ*SQRT( AMAX1( ZERO, ONE - AAPQ*AAPQ ) )
- MXSINJ = AMAX1( MXSINJ, SFMIN )
- END IF
+ CALL SCOPY( M, A( 1, p ), 1,
+ + WORK( N+1 ), 1 )
+ CALL SLASCL( 'G', 0, 0, AAPP, ONE, M,
+ + 1, WORK( N+1 ), LDA,
+ + IERR )
+ CALL SLASCL( 'G', 0, 0, AAQQ, ONE, M,
+ + 1, A( 1, q ), LDA, IERR )
+ TEMP1 = -AAPQ*WORK( p ) / WORK( q )
+ CALL SAXPY( M, TEMP1, WORK( N+1 ), 1,
+ + A( 1, q ), 1 )
+ CALL SLASCL( 'G', 0, 0, ONE, AAQQ, M,
+ + 1, A( 1, q ), LDA, IERR )
+ SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ + ONE-AAPQ*AAPQ ) )
+ MXSINJ = AMAX1( MXSINJ, SFMIN )
+ END IF
* END IF ROTOK THEN ... ELSE
*
* In the case of cancellation in updating SVA(q), SVA(p)
* recompute SVA(q), SVA(p).
*
- IF ( (SVA(q) / AAQQ )**2 .LE. ROOTEPS ) THEN
- IF ((AAQQ .LT. ROOTBIG).AND.(AAQQ .GT. ROOTSFMIN)) THEN
- SVA(q) = SNRM2( M, A(1,q), 1 ) * WORK(q)
- ELSE
- T = ZERO
- AAQQ = ZERO
- CALL SLASSQ( M, A(1,q), 1, T, AAQQ )
- SVA(q) = T * SQRT(AAQQ) * WORK(q)
- END IF
- END IF
- IF ( ( AAPP / AAPP0) .LE. ROOTEPS ) THEN
- IF ((AAPP .LT. ROOTBIG).AND.(AAPP .GT. ROOTSFMIN)) THEN
- AAPP = SNRM2( M, A(1,p), 1 ) * WORK(p)
- ELSE
- T = ZERO
- AAPP = ZERO
- CALL SLASSQ( M, A(1,p), 1, T, AAPP )
- AAPP = T * SQRT(AAPP) * WORK(p)
- END IF
- SVA(p) = AAPP
- END IF
-*
- ELSE
+ IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
+ + THEN
+ IF( ( AAQQ.LT.ROOTBIG ) .AND.
+ + ( AAQQ.GT.ROOTSFMIN ) ) THEN
+ SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
+ + WORK( q )
+ ELSE
+ T = ZERO
+ AAQQ = ZERO
+ CALL SLASSQ( M, A( 1, q ), 1, T,
+ + AAQQ )
+ SVA( q ) = T*SQRT( AAQQ )*WORK( q )
+ END IF
+ END IF
+ IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
+ IF( ( AAPP.LT.ROOTBIG ) .AND.
+ + ( AAPP.GT.ROOTSFMIN ) ) THEN
+ AAPP = SNRM2( M, A( 1, p ), 1 )*
+ + WORK( p )
+ ELSE
+ T = ZERO
+ AAPP = ZERO
+ CALL SLASSQ( M, A( 1, p ), 1, T,
+ + AAPP )
+ AAPP = T*SQRT( AAPP )*WORK( p )
+ END IF
+ SVA( p ) = AAPP
+ END IF
+*
+ ELSE
* A(:,p) and A(:,q) already numerically orthogonal
- IF ( ir1 .EQ. 0 ) NOTROT = NOTROT + 1
+ IF( ir1.EQ.0 )NOTROT = NOTROT + 1
*[RTD] SKIPPED = SKIPPED + 1
- PSKIPPED = PSKIPPED + 1
- END IF
- ELSE
+ PSKIPPED = PSKIPPED + 1
+ END IF
+ ELSE
* A(:,q) is zero column
- IF ( ir1. EQ. 0 ) NOTROT = NOTROT + 1
- PSKIPPED = PSKIPPED + 1
- END IF
+ IF( ir1.EQ.0 )NOTROT = NOTROT + 1
+ PSKIPPED = PSKIPPED + 1
+ END IF
*
- IF ( ( i .LE. SWBAND ) .AND. ( PSKIPPED .GT. ROWSKIP ) ) THEN
- IF ( ir1 .EQ. 0 ) AAPP = - AAPP
- NOTROT = 0
- GO TO 2103
- END IF
+ IF( ( i.LE.SWBAND ) .AND.
+ + ( PSKIPPED.GT.ROWSKIP ) ) THEN
+ IF( ir1.EQ.0 )AAPP = -AAPP
+ NOTROT = 0
+ GO TO 2103
+ END IF
*
- 2002 CONTINUE
+ 2002 CONTINUE
* END q-LOOP
*
- 2103 CONTINUE
+ 2103 CONTINUE
* bailed out of q-loop
*
- SVA(p) = AAPP
+ SVA( p ) = AAPP
*
- ELSE
- SVA(p) = AAPP
- IF ( ( ir1 .EQ. 0 ) .AND. (AAPP .EQ. ZERO) )
- & NOTROT=NOTROT+MIN0(igl+KBL-1,N)-p
- END IF
+ ELSE
+ SVA( p ) = AAPP
+ IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
+ + NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
+ END IF
*
- 2001 CONTINUE
+ 2001 CONTINUE
* end of the p-loop
* end of doing the block ( ibr, ibr )
- 1002 CONTINUE
+ 1002 CONTINUE
* end of ir1-loop
*
* ... go to the off diagonal blocks
*
- igl = ( ibr - 1 ) * KBL + 1
+ igl = ( ibr-1 )*KBL + 1
*
- DO 2010 jbc = ibr + 1, NBL
+ DO 2010 jbc = ibr + 1, NBL
*
- jgl = ( jbc - 1 ) * KBL + 1
+ jgl = ( jbc-1 )*KBL + 1
*
* doing the block at ( ibr, jbc )
*
- IJBLSK = 0
- DO 2100 p = igl, MIN0( igl + KBL - 1, N )
+ IJBLSK = 0
+ DO 2100 p = igl, MIN0( igl+KBL-1, N )
*
- AAPP = SVA(p)
- IF ( AAPP .GT. ZERO ) THEN
+ AAPP = SVA( p )
+ IF( AAPP.GT.ZERO ) THEN
*
- PSKIPPED = 0
+ PSKIPPED = 0
*
- DO 2200 q = jgl, MIN0( jgl + KBL - 1, N )
+ DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
*
- AAQQ = SVA(q)
- IF ( AAQQ .GT. ZERO ) THEN
- AAPP0 = AAPP
+ AAQQ = SVA( q )
+ IF( AAQQ.GT.ZERO ) THEN
+ AAPP0 = AAPP
*
-* -#- M x 2 Jacobi SVD -#-
+* .. M x 2 Jacobi SVD ..
*
* Safe Gram matrix computation
*
- IF ( AAQQ .GE. ONE ) THEN
- IF ( AAPP .GE. AAQQ ) THEN
- ROTOK = ( SMALL*AAPP ) .LE. AAQQ
- ELSE
- ROTOK = ( SMALL*AAQQ ) .LE. AAPP
- END IF
- IF ( AAPP .LT. ( BIG / AAQQ ) ) THEN
- AAPQ = ( SDOT(M, A(1,p), 1, A(1,q), 1 ) *
- & WORK(p) * WORK(q) / AAQQ ) / AAPP
- ELSE
- CALL SCOPY( M, A(1,p), 1, WORK(N+1), 1 )
- CALL SLASCL( 'G', 0, 0, AAPP, WORK(p), M,
- & 1, WORK(N+1), LDA, IERR )
- AAPQ = SDOT( M, WORK(N+1), 1, A(1,q), 1 ) *
- & WORK(q) / AAQQ
- END IF
- ELSE
- IF ( AAPP .GE. AAQQ ) THEN
- ROTOK = AAPP .LE. ( AAQQ / SMALL )
- ELSE
- ROTOK = AAQQ .LE. ( AAPP / SMALL )
- END IF
- IF ( AAPP .GT. ( SMALL / AAQQ ) ) THEN
- AAPQ = ( SDOT( M, A(1,p), 1, A(1,q), 1 ) *
- & WORK(p) * WORK(q) / AAQQ ) / AAPP
- ELSE
- CALL SCOPY( M, A(1,q), 1, WORK(N+1), 1 )
- CALL SLASCL( 'G', 0, 0, AAQQ, WORK(q), M, 1,
- & WORK(N+1), LDA, IERR )
- AAPQ = SDOT(M,WORK(N+1),1,A(1,p),1) * WORK(p) / AAPP
- END IF
- END IF
+ IF( AAQQ.GE.ONE ) THEN
+ IF( AAPP.GE.AAQQ ) THEN
+ ROTOK = ( SMALL*AAPP ).LE.AAQQ
+ ELSE
+ ROTOK = ( SMALL*AAQQ ).LE.AAPP
+ END IF
+ IF( AAPP.LT.( BIG / AAQQ ) ) THEN
+ AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
+ + q ), 1 )*WORK( p )*WORK( q ) /
+ + AAQQ ) / AAPP
+ ELSE
+ CALL SCOPY( M, A( 1, p ), 1,
+ + WORK( N+1 ), 1 )
+ CALL SLASCL( 'G', 0, 0, AAPP,
+ + WORK( p ), M, 1,
+ + WORK( N+1 ), LDA, IERR )
+ AAPQ = SDOT( M, WORK( N+1 ), 1,
+ + A( 1, q ), 1 )*WORK( q ) / AAQQ
+ END IF
+ ELSE
+ IF( AAPP.GE.AAQQ ) THEN
+ ROTOK = AAPP.LE.( AAQQ / SMALL )
+ ELSE
+ ROTOK = AAQQ.LE.( AAPP / SMALL )
+ END IF
+ IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
+ AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
+ + q ), 1 )*WORK( p )*WORK( q ) /
+ + AAQQ ) / AAPP
+ ELSE
+ CALL SCOPY( M, A( 1, q ), 1,
+ + WORK( N+1 ), 1 )
+ CALL SLASCL( 'G', 0, 0, AAQQ,
+ + WORK( q ), M, 1,
+ + WORK( N+1 ), LDA, IERR )
+ AAPQ = SDOT( M, WORK( N+1 ), 1,
+ + A( 1, p ), 1 )*WORK( p ) / AAPP
+ END IF
+ END IF
*
- MXAAPQ = AMAX1( MXAAPQ, ABS(AAPQ) )
+ MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
*
* TO rotate or NOT to rotate, THAT is the question ...
*
- IF ( ABS( AAPQ ) .GT. TOL ) THEN
- NOTROT = 0
+ IF( ABS( AAPQ ).GT.TOL ) THEN
+ NOTROT = 0
*[RTD] ROTATED = ROTATED + 1
- PSKIPPED = 0
- ISWROT = ISWROT + 1
-*
- IF ( ROTOK ) THEN
-*
- AQOAP = AAQQ / AAPP
- APOAQ = AAPP / AAQQ
- THETA = - HALF * ABS( AQOAP - APOAQ ) / AAPQ
- IF ( AAQQ .GT. AAPP0 ) THETA = - THETA
-*
- IF ( ABS( THETA ) .GT. BIGTHETA ) THEN
- T = HALF / THETA
- FASTR(3) = T * WORK(p) / WORK(q)
- FASTR(4) = -T * WORK(q) / WORK(p)
- CALL SROTM( M, A(1,p), 1, A(1,q), 1, FASTR )
- IF ( RSVEC )
- & CALL SROTM( MVL, V(1,p), 1, V(1,q), 1, FASTR )
- SVA(q) = AAQQ*SQRT( AMAX1(ZERO,ONE + T*APOAQ*AAPQ) )
- AAPP = AAPP*SQRT( AMAX1(ZERO,ONE - T*AQOAP*AAPQ) )
- MXSINJ = AMAX1( MXSINJ, ABS(T) )
- ELSE
+ PSKIPPED = 0
+ ISWROT = ISWROT + 1
+*
+ IF( ROTOK ) THEN
+*
+ AQOAP = AAQQ / AAPP
+ APOAQ = AAPP / AAQQ
+ THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
+ IF( AAQQ.GT.AAPP0 )THETA = -THETA
+*
+ IF( ABS( THETA ).GT.BIGTHETA ) THEN
+ T = HALF / THETA
+ FASTR( 3 ) = T*WORK( p ) / WORK( q )
+ FASTR( 4 ) = -T*WORK( q ) /
+ + WORK( p )
+ CALL SROTM( M, A( 1, p ), 1,
+ + A( 1, q ), 1, FASTR )
+ IF( RSVEC )CALL SROTM( MVL,
+ + V( 1, p ), 1,
+ + V( 1, q ), 1,
+ + FASTR )
+ SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ + ONE+T*APOAQ*AAPQ ) )
+ AAPP = AAPP*SQRT( AMAX1( ZERO,
+ + ONE-T*AQOAP*AAPQ ) )
+ MXSINJ = AMAX1( MXSINJ, ABS( T ) )
+ ELSE
*
* .. choose correct signum for THETA and rotate
*
- THSIGN = - SIGN(ONE,AAPQ)
- IF ( AAQQ .GT. AAPP0 ) THSIGN = - THSIGN
- T = ONE / ( THETA + THSIGN*SQRT(ONE+THETA*THETA) )
- CS = SQRT( ONE / ( ONE + T*T ) )
- SN = T * CS
- MXSINJ = AMAX1( MXSINJ, ABS(SN) )
- SVA(q) = AAQQ*SQRT( AMAX1(ZERO, ONE+T*APOAQ*AAPQ) )
- AAPP = AAPP*SQRT( ONE - T*AQOAP*AAPQ)
-*
- APOAQ = WORK(p) / WORK(q)
- AQOAP = WORK(q) / WORK(p)
- IF ( WORK(p) .GE. ONE ) THEN
-*
- IF ( WORK(q) .GE. ONE ) THEN
- FASTR(3) = T * APOAQ
- FASTR(4) = - T * AQOAP
- WORK(p) = WORK(p) * CS
- WORK(q) = WORK(q) * CS
- CALL SROTM( M, A(1,p),1, A(1,q),1, FASTR )
- IF ( RSVEC )
- & CALL SROTM( MVL, V(1,p),1, V(1,q),1, FASTR )
- ELSE
- CALL SAXPY( M, -T*AQOAP, A(1,q),1, A(1,p),1 )
- CALL SAXPY( M, CS*SN*APOAQ, A(1,p),1, A(1,q),1 )
- IF ( RSVEC ) THEN
- CALL SAXPY( MVL, -T*AQOAP, V(1,q),1, V(1,p),1 )
- CALL SAXPY( MVL,CS*SN*APOAQ,V(1,p),1, V(1,q),1 )
- END IF
- WORK(p) = WORK(p) * CS
- WORK(q) = WORK(q) / CS
- END IF
- ELSE
- IF ( WORK(q) .GE. ONE ) THEN
- CALL SAXPY( M, T*APOAQ, A(1,p),1, A(1,q),1 )
- CALL SAXPY( M,-CS*SN*AQOAP, A(1,q),1, A(1,p),1 )
- IF ( RSVEC ) THEN
- CALL SAXPY(MVL,T*APOAQ, V(1,p),1, V(1,q),1 )
- CALL SAXPY(MVL,-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
- CALL SAXPY( M,-T*AQOAP, A(1,q),1,A(1,p),1 )
- CALL SAXPY( M,CS*SN*APOAQ,A(1,p),1,A(1,q),1 )
- WORK(p) = WORK(p) * CS
- WORK(q) = WORK(q) / CS
- IF ( RSVEC ) THEN
- CALL SAXPY( MVL, -T*AQOAP, V(1,q),1,V(1,p),1)
- CALL SAXPY(MVL,CS*SN*APOAQ,V(1,p),1,V(1,q),1)
- END IF
- ELSE
- CALL SAXPY(M, T*APOAQ, A(1,p),1,A(1,q),1)
- CALL SAXPY(M,-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 SAXPY(MVL, T*APOAQ, V(1,p),1,V(1,q),1)
- CALL SAXPY(MVL,-CS*SN*AQOAP,V(1,q),1,V(1,p),1)
- END IF
- END IF
- END IF
- ENDIF
- END IF
-*
- ELSE
- IF ( AAPP .GT. AAQQ ) THEN
- CALL SCOPY( M, A(1,p), 1, WORK(N+1), 1 )
- CALL SLASCL('G',0,0,AAPP,ONE,M,1,WORK(N+1),LDA,IERR)
- CALL SLASCL('G',0,0,AAQQ,ONE,M,1, A(1,q),LDA,IERR)
- TEMP1 = -AAPQ * WORK(p) / WORK(q)
- CALL SAXPY(M,TEMP1,WORK(N+1),1,A(1,q),1)
- CALL SLASCL('G',0,0,ONE,AAQQ,M,1,A(1,q),LDA,IERR)
- SVA(q) = AAQQ*SQRT(AMAX1(ZERO, ONE - AAPQ*AAPQ))
- MXSINJ = AMAX1( MXSINJ, SFMIN )
- ELSE
- CALL SCOPY( M, A(1,q), 1, WORK(N+1), 1 )
- CALL SLASCL('G',0,0,AAQQ,ONE,M,1,WORK(N+1),LDA,IERR)
- CALL SLASCL('G',0,0,AAPP,ONE,M,1, A(1,p),LDA,IERR)
- TEMP1 = -AAPQ * WORK(q) / WORK(p)
- CALL SAXPY(M,TEMP1,WORK(N+1),1,A(1,p),1)
- CALL SLASCL('G',0,0,ONE,AAPP,M,1,A(1,p),LDA,IERR)
- SVA(p) = AAPP*SQRT(AMAX1(ZERO, ONE - AAPQ*AAPQ))
- MXSINJ = AMAX1( MXSINJ, SFMIN )
- END IF
- END IF
+ THSIGN = -SIGN( ONE, AAPQ )
+ IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
+ T = ONE / ( THETA+THSIGN*
+ + SQRT( ONE+THETA*THETA ) )
+ CS = SQRT( ONE / ( ONE+T*T ) )
+ SN = T*CS
+ MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
+ SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ + ONE+T*APOAQ*AAPQ ) )
+ AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
+*
+ APOAQ = WORK( p ) / WORK( q )
+ AQOAP = WORK( q ) / WORK( p )
+ IF( WORK( p ).GE.ONE ) THEN
+*
+ IF( WORK( q ).GE.ONE ) THEN
+ FASTR( 3 ) = T*APOAQ
+ FASTR( 4 ) = -T*AQOAP
+ WORK( p ) = WORK( p )*CS
+ WORK( q ) = WORK( q )*CS
+ CALL SROTM( M, A( 1, p ), 1,
+ + A( 1, q ), 1,
+ + FASTR )
+ IF( RSVEC )CALL SROTM( MVL,
+ + V( 1, p ), 1, V( 1, q ),
+ + 1, FASTR )
+ ELSE
+ CALL SAXPY( M, -T*AQOAP,
+ + A( 1, q ), 1,
+ + A( 1, p ), 1 )
+ CALL SAXPY( M, CS*SN*APOAQ,
+ + A( 1, p ), 1,
+ + A( 1, q ), 1 )
+ IF( RSVEC ) THEN
+ CALL SAXPY( MVL, -T*AQOAP,
+ + V( 1, q ), 1,
+ + V( 1, p ), 1 )
+ CALL SAXPY( MVL,
+ + CS*SN*APOAQ,
+ + V( 1, p ), 1,
+ + V( 1, q ), 1 )
+ END IF
+ WORK( p ) = WORK( p )*CS
+ WORK( q ) = WORK( q ) / CS
+ END IF
+ ELSE
+ IF( WORK( q ).GE.ONE ) THEN
+ CALL SAXPY( M, T*APOAQ,
+ + A( 1, p ), 1,
+ + A( 1, q ), 1 )
+ CALL SAXPY( M, -CS*SN*AQOAP,
+ + A( 1, q ), 1,
+ + A( 1, p ), 1 )
+ IF( RSVEC ) THEN
+ CALL SAXPY( MVL, T*APOAQ,
+ + V( 1, p ), 1,
+ + V( 1, q ), 1 )
+ CALL SAXPY( MVL,
+ + -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
+ CALL SAXPY( M, -T*AQOAP,
+ + A( 1, q ), 1,
+ + A( 1, p ), 1 )
+ CALL SAXPY( M, CS*SN*APOAQ,
+ + A( 1, p ), 1,
+ + A( 1, q ), 1 )
+ WORK( p ) = WORK( p )*CS
+ WORK( q ) = WORK( q ) / CS
+ IF( RSVEC ) THEN
+ CALL SAXPY( MVL,
+ + -T*AQOAP,
+ + V( 1, q ), 1,
+ + V( 1, p ), 1 )
+ CALL SAXPY( MVL,
+ + CS*SN*APOAQ,
+ + V( 1, p ), 1,
+ + V( 1, q ), 1 )
+ END IF
+ ELSE
+ CALL SAXPY( M, T*APOAQ,
+ + A( 1, p ), 1,
+ + A( 1, q ), 1 )
+ CALL SAXPY( M,
+ + -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 SAXPY( MVL,
+ + T*APOAQ, V( 1, p ),
+ + 1, V( 1, q ), 1 )
+ CALL SAXPY( MVL,
+ + -CS*SN*AQOAP,
+ + V( 1, q ), 1,
+ + V( 1, p ), 1 )
+ END IF
+ END IF
+ END IF
+ END IF
+ END IF
+*
+ ELSE
+ IF( AAPP.GT.AAQQ ) THEN
+ CALL SCOPY( M, A( 1, p ), 1,
+ + WORK( N+1 ), 1 )
+ CALL SLASCL( 'G', 0, 0, AAPP, ONE,
+ + M, 1, WORK( N+1 ), LDA,
+ + IERR )
+ CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
+ + M, 1, A( 1, q ), LDA,
+ + IERR )
+ TEMP1 = -AAPQ*WORK( p ) / WORK( q )
+ CALL SAXPY( M, TEMP1, WORK( N+1 ),
+ + 1, A( 1, q ), 1 )
+ CALL SLASCL( 'G', 0, 0, ONE, AAQQ,
+ + M, 1, A( 1, q ), LDA,
+ + IERR )
+ SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+ + ONE-AAPQ*AAPQ ) )
+ MXSINJ = AMAX1( MXSINJ, SFMIN )
+ ELSE
+ CALL SCOPY( M, A( 1, q ), 1,
+ + WORK( N+1 ), 1 )
+ CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
+ + M, 1, WORK( N+1 ), LDA,
+ + IERR )
+ CALL SLASCL( 'G', 0, 0, AAPP, ONE,
+ + M, 1, A( 1, p ), LDA,
+ + IERR )
+ TEMP1 = -AAPQ*WORK( q ) / WORK( p )
+ CALL SAXPY( M, TEMP1, WORK( N+1 ),
+ + 1, A( 1, p ), 1 )
+ CALL SLASCL( 'G', 0, 0, ONE, AAPP,
+ + M, 1, A( 1, p ), LDA,
+ + IERR )
+ SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
+ + ONE-AAPQ*AAPQ ) )
+ MXSINJ = AMAX1( MXSINJ, SFMIN )
+ END IF
+ END IF
* END IF ROTOK THEN ... ELSE
*
* In the case of cancellation in updating SVA(q)
* .. recompute SVA(q)
- IF ( (SVA(q) / AAQQ )**2 .LE. ROOTEPS ) THEN
- IF ((AAQQ .LT. ROOTBIG).AND.(AAQQ .GT. ROOTSFMIN)) THEN
- SVA(q) = SNRM2( M, A(1,q), 1 ) * WORK(q)
- ELSE
- T = ZERO
- AAQQ = ZERO
- CALL SLASSQ( M, A(1,q), 1, T, AAQQ )
- SVA(q) = T * SQRT(AAQQ) * WORK(q)
- END IF
- END IF
- IF ( (AAPP / AAPP0 )**2 .LE. ROOTEPS ) THEN
- IF ((AAPP .LT. ROOTBIG).AND.(AAPP .GT. ROOTSFMIN)) THEN
- AAPP = SNRM2( M, A(1,p), 1 ) * WORK(p)
- ELSE
- T = ZERO
- AAPP = ZERO
- CALL SLASSQ( M, A(1,p), 1, T, AAPP )
- AAPP = T * SQRT(AAPP) * WORK(p)
- END IF
- SVA(p) = AAPP
- END IF
+ IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
+ + THEN
+ IF( ( AAQQ.LT.ROOTBIG ) .AND.
+ + ( AAQQ.GT.ROOTSFMIN ) ) THEN
+ SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
+ + WORK( q )
+ ELSE
+ T = ZERO
+ AAQQ = ZERO
+ CALL SLASSQ( M, A( 1, q ), 1, T,
+ + AAQQ )
+ SVA( q ) = T*SQRT( AAQQ )*WORK( q )
+ END IF
+ END IF
+ IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
+ IF( ( AAPP.LT.ROOTBIG ) .AND.
+ + ( AAPP.GT.ROOTSFMIN ) ) THEN
+ AAPP = SNRM2( M, A( 1, p ), 1 )*
+ + WORK( p )
+ ELSE
+ T = ZERO
+ AAPP = ZERO
+ CALL SLASSQ( M, A( 1, p ), 1, T,
+ + AAPP )
+ AAPP = T*SQRT( AAPP )*WORK( p )
+ END IF
+ SVA( p ) = AAPP
+ END IF
* end of OK rotation
- ELSE
- NOTROT = NOTROT + 1
+ ELSE
+ NOTROT = NOTROT + 1
*[RTD] SKIPPED = SKIPPED + 1
- PSKIPPED = PSKIPPED + 1
- IJBLSK = IJBLSK + 1
- END IF
- ELSE
- NOTROT = NOTROT + 1
- PSKIPPED = PSKIPPED + 1
- IJBLSK = IJBLSK + 1
- END IF
+ PSKIPPED = PSKIPPED + 1
+ IJBLSK = IJBLSK + 1
+ END IF
+ ELSE
+ NOTROT = NOTROT + 1
+ PSKIPPED = PSKIPPED + 1
+ IJBLSK = IJBLSK + 1
+ END IF
*
- IF ( ( i .LE. SWBAND ) .AND. ( IJBLSK .GE. BLSKIP ) ) THEN
- SVA(p) = AAPP
- NOTROT = 0
- GO TO 2011
- END IF
- IF ( ( i .LE. SWBAND ) .AND. ( PSKIPPED .GT. ROWSKIP ) ) THEN
- AAPP = -AAPP
- NOTROT = 0
- GO TO 2203
- END IF
+ IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
+ + THEN
+ SVA( p ) = AAPP
+ NOTROT = 0
+ GO TO 2011
+ END IF
+ IF( ( i.LE.SWBAND ) .AND.
+ + ( PSKIPPED.GT.ROWSKIP ) ) THEN
+ AAPP = -AAPP
+ NOTROT = 0
+ GO TO 2203
+ END IF
*
- 2200 CONTINUE
+ 2200 CONTINUE
* end of the q-loop
- 2203 CONTINUE
+ 2203 CONTINUE
*
- SVA(p) = AAPP
+ SVA( p ) = AAPP
*
- ELSE
+ ELSE
*
- IF ( AAPP .EQ. ZERO ) NOTROT=NOTROT+MIN0(jgl+KBL-1,N)-jgl+1
- IF ( AAPP .LT. ZERO ) NOTROT = 0
+ IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
+ + MIN0( jgl+KBL-1, N ) - jgl + 1
+ IF( AAPP.LT.ZERO )NOTROT = 0
*
- END IF
+ END IF
*
- 2100 CONTINUE
+ 2100 CONTINUE
* end of the p-loop
- 2010 CONTINUE
+ 2010 CONTINUE
* end of the jbc-loop
- 2011 CONTINUE
+ 2011 CONTINUE
*2011 bailed out of the jbc-loop
- DO 2012 p = igl, MIN0( igl + KBL - 1, N )
- SVA(p) = ABS(SVA(p))
- 2012 CONTINUE
+ DO 2012 p = igl, MIN0( igl+KBL-1, N )
+ SVA( p ) = ABS( SVA( p ) )
+ 2012 CONTINUE
***
- 2000 CONTINUE
+ 2000 CONTINUE
*2000 :: end of the ibr-loop
*
* .. update SVA(N)
- IF ((SVA(N) .LT. ROOTBIG).AND.(SVA(N) .GT. ROOTSFMIN)) THEN
- SVA(N) = SNRM2( M, A(1,N), 1 ) * WORK(N)
- ELSE
- T = ZERO
- AAPP = ZERO
- CALL SLASSQ( M, A(1,N), 1, T, AAPP )
- SVA(N) = T * SQRT(AAPP) * WORK(N)
- END IF
+ IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
+ + THEN
+ SVA( N ) = SNRM2( M, A( 1, N ), 1 )*WORK( N )
+ ELSE
+ T = ZERO
+ AAPP = ZERO
+ CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
+ SVA( N ) = T*SQRT( AAPP )*WORK( N )
+ END IF
*
* Additional steering devices
*
- IF ( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
- & ( ISWROT .LE. N ) ) )
- & SWBAND = i
+ IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
+ + ( ISWROT.LE.N ) ) )SWBAND = i
*
- IF ( (i .GT. SWBAND+1) .AND. (MXAAPQ .LT. SQRT(FLOAT(N))*TOL)
- & .AND. (FLOAT(N)*MXAAPQ*MXSINJ .LT. TOL) ) THEN
- GO TO 1994
- END IF
+ IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( FLOAT( N ) )*
+ + TOL ) .AND. ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
+ GO TO 1994
+ END IF
*
- IF ( NOTROT .GE. EMPTSW ) GO TO 1994
+ IF( NOTROT.GE.EMPTSW )GO TO 1994
*
1993 CONTINUE
* end i=1:NSWEEP loop
@@ -1265,80 +1423,81 @@
N2 = 0
N4 = 0
DO 5991 p = 1, N - 1
- q = ISAMAX( N-p+1, SVA(p), 1 ) + p - 1
- IF ( p .NE. q ) THEN
- TEMP1 = SVA(p)
- SVA(p) = SVA(q)
- SVA(q) = TEMP1
- TEMP1 = WORK(p)
- WORK(p) = WORK(q)
- WORK(q) = TEMP1
- CALL SSWAP( M, A(1,p), 1, A(1,q), 1 )
- IF ( RSVEC ) CALL SSWAP( MVL, V(1,p), 1, V(1,q), 1 )
+ q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
+ IF( p.NE.q ) THEN
+ TEMP1 = SVA( p )
+ SVA( p ) = SVA( q )
+ SVA( q ) = TEMP1
+ TEMP1 = WORK( p )
+ WORK( p ) = WORK( q )
+ WORK( q ) = TEMP1
+ CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
+ IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
END IF
- IF ( SVA(p) .NE. ZERO ) THEN
+ IF( SVA( p ).NE.ZERO ) THEN
N4 = N4 + 1
- IF ( SVA(p)*SCALE .GT. SFMIN ) N2 = N2 + 1
+ IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1
END IF
5991 CONTINUE
- IF ( SVA(N) .NE. ZERO ) THEN
+ IF( SVA( N ).NE.ZERO ) THEN
N4 = N4 + 1
- IF ( SVA(N)*SCALE .GT. SFMIN ) N2 = N2 + 1
+ IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1
END IF
*
* Normalize the left singular vectors.
*
- IF ( LSVEC .OR. UCTOL ) THEN
+ IF( LSVEC .OR. UCTOL ) THEN
DO 1998 p = 1, N2
- CALL SSCAL( M, WORK(p) / SVA(p), A(1,p), 1 )
+ CALL SSCAL( M, WORK( p ) / SVA( p ), A( 1, p ), 1 )
1998 CONTINUE
END IF
*
* Scale the product of Jacobi rotations (assemble the fast rotations).
*
- IF ( RSVEC ) THEN
- IF ( APPLV ) THEN
+ IF( RSVEC ) THEN
+ IF( APPLV ) THEN
DO 2398 p = 1, N
- CALL SSCAL( MVL, WORK(p), V(1,p), 1 )
+ CALL SSCAL( MVL, WORK( p ), V( 1, p ), 1 )
2398 CONTINUE
ELSE
DO 2399 p = 1, N
- TEMP1 = ONE / SNRM2(MVL, V(1,p), 1 )
- CALL SSCAL( MVL, TEMP1, V(1,p), 1 )
+ TEMP1 = ONE / SNRM2( MVL, V( 1, p ), 1 )
+ CALL SSCAL( MVL, TEMP1, V( 1, p ), 1 )
2399 CONTINUE
END IF
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( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
+ + SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.
+ + ( SFMIN / SCALE ) ) ) ) THEN
DO 2400 p = 1, N
- SVA(p) = SCALE*SVA(p)
+ SVA( p ) = SCALE*SVA( p )
2400 CONTINUE
SCALE = ONE
END IF
*
- WORK(1) = SCALE
+ WORK( 1 ) = SCALE
* The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE
* then some of the singular values may overflow or underflow and
* the spectrum is given in this factored representation.
*
- WORK(2) = FLOAT(N4)
+ WORK( 2 ) = FLOAT( N4 )
* N4 is the number of computed nonzero singular values of A.
*
- WORK(3) = FLOAT(N2)
+ WORK( 3 ) = FLOAT( N2 )
* N2 is the number of singular values of A greater than SFMIN.
* If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
* that may carry some information.
*
- WORK(4) = FLOAT(i)
+ WORK( 4 ) = FLOAT( i )
* i is the index of the last sweep before declaring convergence.
*
- WORK(5) = MXAAPQ
+ WORK( 5 ) = MXAAPQ
* MXAAPQ is the largest absolute value of scaled pivots in the
* last sweep
*
- WORK(6) = MXSINJ
+ WORK( 6 ) = MXSINJ
* MXSINJ is the largest absolute value of the sines of Jacobi angles
* in the last sweep
*
@@ -1347,4 +1506,3 @@
* .. END OF SGESVJ
* ..
END
-*
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-21 17:41:55 +0000