Annotation of rpl/lapack/lapack/zlargv.f, revision 1.4
1.1 bertrand 1: SUBROUTINE ZLARGV( N, X, INCX, Y, INCY, C, INCC )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER INCC, INCX, INCY, N
10: * ..
11: * .. Array Arguments ..
12: DOUBLE PRECISION C( * )
13: COMPLEX*16 X( * ), Y( * )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * ZLARGV generates a vector of complex plane rotations with real
20: * cosines, determined by elements of the complex vectors x and y.
21: * For i = 1,2,...,n
22: *
23: * ( c(i) s(i) ) ( x(i) ) = ( r(i) )
24: * ( -conjg(s(i)) c(i) ) ( y(i) ) = ( 0 )
25: *
26: * where c(i)**2 + ABS(s(i))**2 = 1
27: *
28: * The following conventions are used (these are the same as in ZLARTG,
29: * but differ from the BLAS1 routine ZROTG):
30: * If y(i)=0, then c(i)=1 and s(i)=0.
31: * If x(i)=0, then c(i)=0 and s(i) is chosen so that r(i) is real.
32: *
33: * Arguments
34: * =========
35: *
36: * N (input) INTEGER
37: * The number of plane rotations to be generated.
38: *
39: * X (input/output) COMPLEX*16 array, dimension (1+(N-1)*INCX)
40: * On entry, the vector x.
41: * On exit, x(i) is overwritten by r(i), for i = 1,...,n.
42: *
43: * INCX (input) INTEGER
44: * The increment between elements of X. INCX > 0.
45: *
46: * Y (input/output) COMPLEX*16 array, dimension (1+(N-1)*INCY)
47: * On entry, the vector y.
48: * On exit, the sines of the plane rotations.
49: *
50: * INCY (input) INTEGER
51: * The increment between elements of Y. INCY > 0.
52: *
53: * C (output) DOUBLE PRECISION array, dimension (1+(N-1)*INCC)
54: * The cosines of the plane rotations.
55: *
56: * INCC (input) INTEGER
57: * The increment between elements of C. INCC > 0.
58: *
59: * Further Details
60: * ======= =======
61: *
62: * 6-6-96 - Modified with a new algorithm by W. Kahan and J. Demmel
63: *
64: * This version has a few statements commented out for thread safety
65: * (machine parameters are computed on each entry). 10 feb 03, SJH.
66: *
67: * =====================================================================
68: *
69: * .. Parameters ..
70: DOUBLE PRECISION TWO, ONE, ZERO
71: PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
72: COMPLEX*16 CZERO
73: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
74: * ..
75: * .. Local Scalars ..
76: * LOGICAL FIRST
77:
78: INTEGER COUNT, I, IC, IX, IY, J
79: DOUBLE PRECISION CS, D, DI, DR, EPS, F2, F2S, G2, G2S, SAFMIN,
80: $ SAFMN2, SAFMX2, SCALE
81: COMPLEX*16 F, FF, FS, G, GS, R, SN
82: * ..
83: * .. External Functions ..
84: DOUBLE PRECISION DLAMCH, DLAPY2
85: EXTERNAL DLAMCH, DLAPY2
86: * ..
87: * .. Intrinsic Functions ..
88: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, LOG,
89: $ MAX, SQRT
90: * ..
91: * .. Statement Functions ..
92: DOUBLE PRECISION ABS1, ABSSQ
93: * ..
94: * .. Save statement ..
95: * SAVE FIRST, SAFMX2, SAFMIN, SAFMN2
96: * ..
97: * .. Data statements ..
98: * DATA FIRST / .TRUE. /
99: * ..
100: * .. Statement Function definitions ..
101: ABS1( FF ) = MAX( ABS( DBLE( FF ) ), ABS( DIMAG( FF ) ) )
102: ABSSQ( FF ) = DBLE( FF )**2 + DIMAG( FF )**2
103: * ..
104: * .. Executable Statements ..
105: *
106: * IF( FIRST ) THEN
107: * FIRST = .FALSE.
108: SAFMIN = DLAMCH( 'S' )
109: EPS = DLAMCH( 'E' )
110: SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
111: $ LOG( DLAMCH( 'B' ) ) / TWO )
112: SAFMX2 = ONE / SAFMN2
113: * END IF
114: IX = 1
115: IY = 1
116: IC = 1
117: DO 60 I = 1, N
118: F = X( IX )
119: G = Y( IY )
120: *
121: * Use identical algorithm as in ZLARTG
122: *
123: SCALE = MAX( ABS1( F ), ABS1( G ) )
124: FS = F
125: GS = G
126: COUNT = 0
127: IF( SCALE.GE.SAFMX2 ) THEN
128: 10 CONTINUE
129: COUNT = COUNT + 1
130: FS = FS*SAFMN2
131: GS = GS*SAFMN2
132: SCALE = SCALE*SAFMN2
133: IF( SCALE.GE.SAFMX2 )
134: $ GO TO 10
135: ELSE IF( SCALE.LE.SAFMN2 ) THEN
136: IF( G.EQ.CZERO ) THEN
137: CS = ONE
138: SN = CZERO
139: R = F
140: GO TO 50
141: END IF
142: 20 CONTINUE
143: COUNT = COUNT - 1
144: FS = FS*SAFMX2
145: GS = GS*SAFMX2
146: SCALE = SCALE*SAFMX2
147: IF( SCALE.LE.SAFMN2 )
148: $ GO TO 20
149: END IF
150: F2 = ABSSQ( FS )
151: G2 = ABSSQ( GS )
152: IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
153: *
154: * This is a rare case: F is very small.
155: *
156: IF( F.EQ.CZERO ) THEN
157: CS = ZERO
158: R = DLAPY2( DBLE( G ), DIMAG( G ) )
159: * Do complex/real division explicitly with two real
160: * divisions
161: D = DLAPY2( DBLE( GS ), DIMAG( GS ) )
162: SN = DCMPLX( DBLE( GS ) / D, -DIMAG( GS ) / D )
163: GO TO 50
164: END IF
165: F2S = DLAPY2( DBLE( FS ), DIMAG( FS ) )
166: * G2 and G2S are accurate
167: * G2 is at least SAFMIN, and G2S is at least SAFMN2
168: G2S = SQRT( G2 )
169: * Error in CS from underflow in F2S is at most
170: * UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
171: * If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
172: * and so CS .lt. sqrt(SAFMIN)
173: * If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
174: * and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
175: * Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
176: CS = F2S / G2S
177: * Make sure abs(FF) = 1
178: * Do complex/real division explicitly with 2 real divisions
179: IF( ABS1( F ).GT.ONE ) THEN
180: D = DLAPY2( DBLE( F ), DIMAG( F ) )
181: FF = DCMPLX( DBLE( F ) / D, DIMAG( F ) / D )
182: ELSE
183: DR = SAFMX2*DBLE( F )
184: DI = SAFMX2*DIMAG( F )
185: D = DLAPY2( DR, DI )
186: FF = DCMPLX( DR / D, DI / D )
187: END IF
188: SN = FF*DCMPLX( DBLE( GS ) / G2S, -DIMAG( GS ) / G2S )
189: R = CS*F + SN*G
190: ELSE
191: *
192: * This is the most common case.
193: * Neither F2 nor F2/G2 are less than SAFMIN
194: * F2S cannot overflow, and it is accurate
195: *
196: F2S = SQRT( ONE+G2 / F2 )
197: * Do the F2S(real)*FS(complex) multiply with two real
198: * multiplies
199: R = DCMPLX( F2S*DBLE( FS ), F2S*DIMAG( FS ) )
200: CS = ONE / F2S
201: D = F2 + G2
202: * Do complex/real division explicitly with two real divisions
203: SN = DCMPLX( DBLE( R ) / D, DIMAG( R ) / D )
204: SN = SN*DCONJG( GS )
205: IF( COUNT.NE.0 ) THEN
206: IF( COUNT.GT.0 ) THEN
207: DO 30 J = 1, COUNT
208: R = R*SAFMX2
209: 30 CONTINUE
210: ELSE
211: DO 40 J = 1, -COUNT
212: R = R*SAFMN2
213: 40 CONTINUE
214: END IF
215: END IF
216: END IF
217: 50 CONTINUE
218: C( IC ) = CS
219: Y( IY ) = SN
220: X( IX ) = R
221: IC = IC + INCC
222: IY = IY + INCY
223: IX = IX + INCX
224: 60 CONTINUE
225: RETURN
226: *
227: * End of ZLARGV
228: *
229: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlargv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack