diff options
| author | julie <julielangou@users.noreply.github.com> | 2010-08-25 15:59:11 +0000 |
|---|---|---|
| committer | julie <julielangou@users.noreply.github.com> | 2010-08-25 15:59:11 +0000 |
| commit | f327dddb7f08d3738dc267f46b6fb06849b0a6a4 (patch) | |
| tree | bcc22a22ce63312e34046f41084825b545c1e0f3 /SRC/dgsvj0.f | |
| parent | 150cd6266c7d887da544f7f23bfd7fd610044f66 (diff) | |
| download | lapack-f327dddb7f08d3738dc267f46b6fb06849b0a6a4.tar.gz lapack-f327dddb7f08d3738dc267f46b6fb06849b0a6a4.tar.bz2 lapack-f327dddb7f08d3738dc267f46b6fb06849b0a6a4.zip | |
Included bug fix provided by Zlatco on Jacobi SVD
Email from Zlatco on August 24th 2010:
The problem that was reported (with zero matrix) is caused by bad initialization to xLASSQ.
It should be ZERO, ONE and not ZERO, ZERO. In fact, I had it ZERO, ONE throughout
the complete development of the code and decided to change it to ZERO, ZERO a the very
end to make it "more elegant". That was stupid, because xLASSQ does not touch those
variables in case of zero vector, leaving scaling at ZERO, and in the nonzero case the scaling
is between ONE and SQRT(N). So, in case of zero vector, division by a variable that
is normally bigger than ONE causes division by zero.
I have corrected that and few other things, stress tested the code and it should be OK now.
README:
i) In xgejsv.f and xgesvj.f input parameters SCALE and
SUMSQ in xlassq.f are now initially set as SCALE = ZERO, SUMSQ=ONE.
Setting them both to zero (without carefully reading xlassq.f) caused
problems with exactly zero columns.
ii) There was a problem in the branch that computes only SIGMA and U of a
rank deficient matrix. The computed numerical rank (NR) was incorrectly
written as N in parameter lists of the corresponding calls.
iii) In xgsvj0.f, xgsvj1.f testing the input parameters is changed to prevent
unnecessarily negative INFO in some situations.
iv) Minor changes, renaming some variables etc.
Diffstat (limited to 'SRC/dgsvj0.f')
| -rw-r--r-- | SRC/dgsvj0.f | 25 |
1 files changed, 13 insertions, 12 deletions
diff --git a/SRC/dgsvj0.f b/SRC/dgsvj0.f index 2df0181c..0312690a 100644 --- a/SRC/dgsvj0.f +++ b/SRC/dgsvj0.f @@ -183,9 +183,10 @@ INFO = -3 ELSE IF( LDA.LT.M ) THEN INFO = -5 - ELSE IF( MV.LT.0 ) THEN + ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN INFO = -8 - ELSE IF( LDV.LT.M ) THEN + ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. + & ( APPLV.AND.( LDV.LT.MV ) ) ) THEN INFO = -10 ELSE IF( TOL.LE.EPS ) THEN INFO = -13 @@ -307,7 +308,7 @@ SVA( p ) = DNRM2( M, A( 1, p ), 1 )*D( p ) ELSE TEMP1 = ZERO - AAPP = ZERO + AAPP = ONE CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) SVA( p ) = TEMP1*DSQRT( AAPP )*D( p ) END IF @@ -391,8 +392,8 @@ + FASTR ) SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, + ONE+T*APOAQ*AAPQ ) ) - AAPP = AAPP*DSQRT( ONE-T*AQOAP* - + AAPQ ) + AAPP = AAPP*DSQRT( DMAX1( ZERO, + + ONE-T*AQOAP*AAPQ ) ) MXSINJ = DMAX1( MXSINJ, DABS( T ) ) * ELSE @@ -535,7 +536,7 @@ + D( q ) ELSE T = ZERO - AAQQ = ZERO + AAQQ = ONE CALL DLASSQ( M, A( 1, q ), 1, T, + AAQQ ) SVA( q ) = T*DSQRT( AAQQ )*D( q ) @@ -548,7 +549,7 @@ + D( p ) ELSE T = ZERO - AAPP = ZERO + AAPP = ONE CALL DLASSQ( M, A( 1, p ), 1, T, + AAPP ) AAPP = T*DSQRT( AAPP )*D( p ) @@ -707,8 +708,8 @@ MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, + ONE+T*APOAQ*AAPQ ) ) - AAPP = AAPP*DSQRT( ONE-T*AQOAP* - + AAPQ ) + AAPP = AAPP*DSQRT( DMAX1( ZERO, + + ONE-T*AQOAP*AAPQ ) ) * APOAQ = D( p ) / D( q ) AQOAP = D( q ) / D( p ) @@ -856,7 +857,7 @@ + D( q ) ELSE T = ZERO - AAQQ = ZERO + AAQQ = ONE CALL DLASSQ( M, A( 1, q ), 1, T, + AAQQ ) SVA( q ) = T*DSQRT( AAQQ )*D( q ) @@ -869,7 +870,7 @@ + D( p ) ELSE T = ZERO - AAPP = ZERO + AAPP = ONE CALL DLASSQ( M, A( 1, p ), 1, T, + AAPP ) AAPP = T*DSQRT( AAPP )*D( p ) @@ -932,7 +933,7 @@ SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N ) ELSE T = ZERO - AAPP = ZERO + AAPP = ONE CALL DLASSQ( M, A( 1, N ), 1, T, AAPP ) SVA( N ) = T*DSQRT( AAPP )*D( N ) END IF |