Mercurial > octave-nkf

annotate libcruft/lapack/dlartg.f @ 2329:30c606bec7a8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-07-19 01:29:05 by jwe] Initial revision
author jwe
date Fri, 19 Jul 1996 01:29:55 +0000
parents
children 15cddaacbc2d
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
1 SUBROUTINE DLARTG( F, G, CS, SN, R )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
2 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
3 * -- LAPACK auxiliary routine (version 2.0) --
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
5 * Courant Institute, Argonne National Lab, and Rice University
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
6 * September 30, 1994
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
7 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
8 * .. Scalar Arguments ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
9 DOUBLE PRECISION CS, F, G, R, SN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
10 * ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
11 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
12 * Purpose
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
13 * =======
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
14 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
15 * DLARTG generate a plane rotation so that
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
16 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
17 * [ CS SN ] . [ F ] = [ R ] where CS**2 + SN**2 = 1.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
18 * [ -SN CS ] [ G ] [ 0 ]
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
19 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
20 * This is a slower, more accurate version of the BLAS1 routine DROTG,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
21 * with the following other differences:
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
22 * F and G are unchanged on return.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
23 * If G=0, then CS=1 and SN=0.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
24 * If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
25 * floating point operations (saves work in DBDSQR when
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
26 * there are zeros on the diagonal).
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
27 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
28 * If F exceeds G in magnitude, CS will be positive.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
29 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
30 * Arguments
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
31 * =========
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
32 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
33 * F (input) DOUBLE PRECISION
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
34 * The first component of vector to be rotated.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
35 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
36 * G (input) DOUBLE PRECISION
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
37 * The second component of vector to be rotated.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
38 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
39 * CS (output) DOUBLE PRECISION
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
40 * The cosine of the rotation.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
41 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
42 * SN (output) DOUBLE PRECISION
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
43 * The sine of the rotation.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
44 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
45 * R (output) DOUBLE PRECISION
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
46 * The nonzero component of the rotated vector.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
47 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
48 * =====================================================================
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
49 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
50 * .. Parameters ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
51 DOUBLE PRECISION ZERO
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
52 PARAMETER ( ZERO = 0.0D0 )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
53 DOUBLE PRECISION ONE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
54 PARAMETER ( ONE = 1.0D0 )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
55 DOUBLE PRECISION TWO
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
56 PARAMETER ( TWO = 2.0D0 )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
57 * ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
58 * .. Local Scalars ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
59 LOGICAL FIRST
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
60 INTEGER COUNT, I
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
61 DOUBLE PRECISION EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
62 * ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
63 * .. External Functions ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
64 DOUBLE PRECISION DLAMCH
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
65 EXTERNAL DLAMCH
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
66 * ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
67 * .. Intrinsic Functions ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
68 INTRINSIC ABS, INT, LOG, MAX, SQRT
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
69 * ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
70 * .. Save statement ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
71 SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
72 * ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
73 * .. Data statements ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
74 DATA FIRST / .TRUE. /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
75 * ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
76 * .. Executable Statements ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
77 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
78 IF( FIRST ) THEN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
79 FIRST = .FALSE.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
80 SAFMIN = DLAMCH( 'S' )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
81 EPS = DLAMCH( 'E' )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
82 SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
83 $ LOG( DLAMCH( 'B' ) ) / TWO )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
84 SAFMX2 = ONE / SAFMN2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
85 END IF
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
86 IF( G.EQ.ZERO ) THEN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
87 CS = ONE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
88 SN = ZERO
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
89 R = F
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
90 ELSE IF( F.EQ.ZERO ) THEN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
91 CS = ZERO
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
92 SN = ONE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
93 R = G
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
94 ELSE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
95 F1 = F
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
96 G1 = G
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
97 SCALE = MAX( ABS( F1 ), ABS( G1 ) )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
98 IF( SCALE.GE.SAFMX2 ) THEN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
99 COUNT = 0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
100 10 CONTINUE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
101 COUNT = COUNT + 1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
102 F1 = F1*SAFMN2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
103 G1 = G1*SAFMN2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
104 SCALE = MAX( ABS( F1 ), ABS( G1 ) )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
105 IF( SCALE.GE.SAFMX2 )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
106 $ GO TO 10
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
107 R = SQRT( F1**2+G1**2 )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
108 CS = F1 / R
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
109 SN = G1 / R
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
110 DO 20 I = 1, COUNT
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
111 R = R*SAFMX2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
112 20 CONTINUE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
113 ELSE IF( SCALE.LE.SAFMN2 ) THEN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
114 COUNT = 0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
115 30 CONTINUE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
116 COUNT = COUNT + 1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
117 F1 = F1*SAFMX2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
118 G1 = G1*SAFMX2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
119 SCALE = MAX( ABS( F1 ), ABS( G1 ) )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
120 IF( SCALE.LE.SAFMN2 )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
121 $ GO TO 30
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
122 R = SQRT( F1**2+G1**2 )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
123 CS = F1 / R
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
124 SN = G1 / R
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
125 DO 40 I = 1, COUNT
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
126 R = R*SAFMN2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
127 40 CONTINUE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
128 ELSE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
129 R = SQRT( F1**2+G1**2 )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
130 CS = F1 / R
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
131 SN = G1 / R
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
132 END IF
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
133 IF( ABS( F ).GT.ABS( G ) .AND. CS.LT.ZERO ) THEN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
134 CS = -CS
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
135 SN = -SN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
136 R = -R
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
137 END IF
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
138 END IF
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
139 RETURN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
140 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
141 * End of DLARTG
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
142 *
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
143 END