diff options
| author | brian <brian@8a072113-8704-0410-8d35-dd094bca7971> | 2010-11-03 23:02:29 +0000 |
|---|---|---|
| committer | brian <brian@8a072113-8704-0410-8d35-dd094bca7971> | 2010-11-03 23:02:29 +0000 |
| commit | 4ca2feaf79883558f849f792f6813819da97a821 (patch) | |
| tree | 7079f3949a0356cd2914ab4984e928ef2ebf1b8e /SRC/slartgp.f | |
| parent | 1237a0d5b7f033a117062f78bf055026928af9ec (diff) | |
| download | lapack-4ca2feaf79883558f849f792f6813819da97a821.tar.gz lapack-4ca2feaf79883558f849f792f6813819da97a821.tar.bz2 lapack-4ca2feaf79883558f849f792f6813819da97a821.zip | |
Added CS decomposition source files to SRC/
Diffstat (limited to 'SRC/slartgp.f')
| -rw-r--r-- | SRC/slartgp.f | 148 |
1 files changed, 148 insertions, 0 deletions
diff --git a/SRC/slartgp.f b/SRC/slartgp.f new file mode 100644 index 00000000..9bb635e3 --- /dev/null +++ b/SRC/slartgp.f @@ -0,0 +1,148 @@ + SUBROUTINE SLARTGP( F, G, CS, SN, R ) +* +* Originally SLARTG +* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2006 +* +* Adapted to SLARTGP by Brian Sutton +* July 2010 +* +* .. Scalar Arguments .. + REAL CS, F, G, R, SN +* .. +* +* Purpose +* ======= +* +* SLARTGP generates a plane rotation so that +* +* [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1. +* [ -SN CS ] [ G ] [ 0 ] +* +* This is a slower, more accurate version of the Level 1 BLAS routine SROTG, +* with the following other differences: +* F and G are unchanged on return. +* If G=0, then CS=(+/-)1 and SN=0. +* If F=0 and (G .ne. 0), then CS=0 and SN=(+/-)1. +* +* The sign is chosen so that R >= 0. +* +* Arguments +* ========= +* +* F (input) REAL +* The first component of vector to be rotated. +* +* G (input) REAL +* The second component of vector to be rotated. +* +* CS (output) REAL +* The cosine of the rotation. +* +* SN (output) REAL +* The sine of the rotation. +* +* R (output) REAL +* The nonzero component of the rotated vector. +* +* This version has a few statements commented out for thread safety +* (machine parameters are computed on each entry). 10 feb 03, SJH. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0E0 ) + REAL ONE + PARAMETER ( ONE = 1.0E0 ) + REAL TWO + PARAMETER ( TWO = 2.0E0 ) +* .. +* .. Local Scalars .. +* LOGICAL FIRST + INTEGER COUNT, I + REAL EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE +* .. +* .. External Functions .. + REAL SLAMCH + EXTERNAL SLAMCH +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, INT, LOG, MAX, SIGN, SQRT +* .. +* .. Save statement .. +* SAVE FIRST, SAFMX2, SAFMIN, SAFMN2 +* .. +* .. Data statements .. +* DATA FIRST / .TRUE. / +* .. +* .. Executable Statements .. +* +* IF( FIRST ) THEN + SAFMIN = SLAMCH( 'S' ) + EPS = SLAMCH( 'E' ) + SAFMN2 = SLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / + $ LOG( SLAMCH( 'B' ) ) / TWO ) + SAFMX2 = ONE / SAFMN2 +* FIRST = .FALSE. +* END IF + IF( G.EQ.ZERO ) THEN + CS = SIGN( ONE, F ) + SN = ZERO + R = ABS( F ) + ELSE IF( F.EQ.ZERO ) THEN + CS = ZERO + SN = SIGN( ONE, G ) + R = ABS( G ) + ELSE + F1 = F + G1 = G + SCALE = MAX( ABS( F1 ), ABS( G1 ) ) + IF( SCALE.GE.SAFMX2 ) THEN + COUNT = 0 + 10 CONTINUE + COUNT = COUNT + 1 + F1 = F1*SAFMN2 + G1 = G1*SAFMN2 + SCALE = MAX( ABS( F1 ), ABS( G1 ) ) + IF( SCALE.GE.SAFMX2 ) + $ GO TO 10 + R = SQRT( F1**2+G1**2 ) + CS = F1 / R + SN = G1 / R + DO 20 I = 1, COUNT + R = R*SAFMX2 + 20 CONTINUE + ELSE IF( SCALE.LE.SAFMN2 ) THEN + COUNT = 0 + 30 CONTINUE + COUNT = COUNT + 1 + F1 = F1*SAFMX2 + G1 = G1*SAFMX2 + SCALE = MAX( ABS( F1 ), ABS( G1 ) ) + IF( SCALE.LE.SAFMN2 ) + $ GO TO 30 + R = SQRT( F1**2+G1**2 ) + CS = F1 / R + SN = G1 / R + DO 40 I = 1, COUNT + R = R*SAFMN2 + 40 CONTINUE + ELSE + R = SQRT( F1**2+G1**2 ) + CS = F1 / R + SN = G1 / R + END IF + IF( R.LT.ZERO ) THEN + CS = -CS + SN = -SN + R = -R + END IF + END IF + RETURN +* +* End of SLARTG +* + END |