source: trunk/SOURCES/BLAS/slaic1.f @ 22

Last change on this file since 22 was 22, checked in by roche, 9 years ago

Petites adaptations diverses du code pour compilation en gfortran. Ajout d un Makefile flexible a option pour choisir ifort ou gfortran.
File size: 9.8 KB
Line 
1*> \brief \b SLAIC1
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at 
6*            http://www.netlib.org/lapack/explore-html/ 
7*
8*> \htmlonly
9*> Download SLAIC1 + dependencies 
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaic1.f"> 
11*> [TGZ]</a> 
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaic1.f"> 
13*> [ZIP]</a> 
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaic1.f"> 
15*> [TXT]</a>
16*> \endhtmlonly 
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
22* 
23*       .. Scalar Arguments ..
24*       INTEGER            J, JOB
25*       REAL               C, GAMMA, S, SEST, SESTPR
26*       ..
27*       .. Array Arguments ..
28*       REAL               W( J ), X( J )
29*       ..
30* 
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> SLAIC1 applies one step of incremental condition estimation in
38*> its simplest version:
39*>
40*> Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
41*> lower triangular matrix L, such that
42*>          twonorm(L*x) = sest
43*> Then SLAIC1 computes sestpr, s, c such that
44*> the vector
45*>                 [ s*x ]
46*>          xhat = [  c  ]
47*> is an approximate singular vector of
48*>                 [ L      0  ]
49*>          Lhat = [ w**T gamma ]
50*> in the sense that
51*>          twonorm(Lhat*xhat) = sestpr.
52*>
53*> Depending on JOB, an estimate for the largest or smallest singular
54*> value is computed.
55*>
56*> Note that [s c]**T and sestpr**2 is an eigenpair of the system
57*>
58*>     diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
59*>                                           [ gamma ]
60*>
61*> where  alpha =  x**T*w.
62*> \endverbatim
63*
64*  Arguments:
65*  ==========
66*
67*> \param[in] JOB
68*> \verbatim
69*>          JOB is INTEGER
70*>          = 1: an estimate for the largest singular value is computed.
71*>          = 2: an estimate for the smallest singular value is computed.
72*> \endverbatim
73*>
74*> \param[in] J
75*> \verbatim
76*>          J is INTEGER
77*>          Length of X and W
78*> \endverbatim
79*>
80*> \param[in] X
81*> \verbatim
82*>          X is REAL array, dimension (J)
83*>          The j-vector x.
84*> \endverbatim
85*>
86*> \param[in] SEST
87*> \verbatim
88*>          SEST is REAL
89*>          Estimated singular value of j by j matrix L
90*> \endverbatim
91*>
92*> \param[in] W
93*> \verbatim
94*>          W is REAL array, dimension (J)
95*>          The j-vector w.
96*> \endverbatim
97*>
98*> \param[in] GAMMA
99*> \verbatim
100*>          GAMMA is REAL
101*>          The diagonal element gamma.
102*> \endverbatim
103*>
104*> \param[out] SESTPR
105*> \verbatim
106*>          SESTPR is REAL
107*>          Estimated singular value of (j+1) by (j+1) matrix Lhat.
108*> \endverbatim
109*>
110*> \param[out] S
111*> \verbatim
112*>          S is REAL
113*>          Sine needed in forming xhat.
114*> \endverbatim
115*>
116*> \param[out] C
117*> \verbatim
118*>          C is REAL
119*>          Cosine needed in forming xhat.
120*> \endverbatim
121*
122*  Authors:
123*  ========
124*
125*> \author Univ. of Tennessee 
126*> \author Univ. of California Berkeley 
127*> \author Univ. of Colorado Denver 
128*> \author NAG Ltd. 
129*
130*> \date November 2011
131*
132*> \ingroup realOTHERauxiliary
133*
134*  =====================================================================
135      SUBROUTINE SLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C )
136*
137*  -- LAPACK auxiliary routine (version 3.4.0) --
138*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
139*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140*     November 2011
141*
142*     .. Scalar Arguments ..
143      INTEGER            J, JOB
144      REAL               C, GAMMA, S, SEST, SESTPR
145*     ..
146*     .. Array Arguments ..
147      REAL               W( J ), X( J )
148*     ..
149*
150*  =====================================================================
151*
152*     .. Parameters ..
153      REAL               ZERO, ONE, TWO
154      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
155      REAL               HALF, FOUR
156      PARAMETER          ( HALF = 0.5E0, FOUR = 4.0E0 )
157*     ..
158*     .. Local Scalars ..
159      REAL               ABSALP, ABSEST, ABSGAM, ALPHA, B, COSINE, EPS,
160     $                   NORMA, S1, S2, SINE, T, TEST, TMP, ZETA1, ZETA2
161*     ..
162*     .. Intrinsic Functions ..
163      INTRINSIC          ABS, MAX, SIGN, SQRT
164*     ..
165*     .. External Functions ..
166      REAL               SDOT, SLAMCH
167      EXTERNAL           SDOT, SLAMCH
168*     ..
169*     .. Executable Statements ..
170*
171      EPS = SLAMCH( 'Epsilon' )
172      ALPHA = SDOT( J, X, 1, W, 1 )
173*
174      ABSALP = ABS( ALPHA )
175      ABSGAM = ABS( GAMMA )
176      ABSEST = ABS( SEST )
177*
178      IF( JOB.EQ.1 ) THEN
179*
180*        Estimating largest singular value
181*
182*        special cases
183*
184         IF( SEST.EQ.ZERO ) THEN
185            S1 = MAX( ABSGAM, ABSALP )
186            IF( S1.EQ.ZERO ) THEN
187               S = ZERO
188               C = ONE
189               SESTPR = ZERO
190            ELSE
191               S = ALPHA / S1
192               C = GAMMA / S1
193               TMP = SQRT( S*S+C*C )
194               S = S / TMP
195               C = C / TMP
196               SESTPR = S1*TMP
197            END IF
198            RETURN
199         ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
200            S = ONE
201            C = ZERO
202            TMP = MAX( ABSEST, ABSALP )
203            S1 = ABSEST / TMP
204            S2 = ABSALP / TMP
205            SESTPR = TMP*SQRT( S1*S1+S2*S2 )
206            RETURN
207         ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
208            S1 = ABSGAM
209            S2 = ABSEST
210            IF( S1.LE.S2 ) THEN
211               S = ONE
212               C = ZERO
213               SESTPR = S2
214            ELSE
215               S = ZERO
216               C = ONE
217               SESTPR = S1
218            END IF
219            RETURN
220         ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
221            S1 = ABSGAM
222            S2 = ABSALP
223            IF( S1.LE.S2 ) THEN
224               TMP = S1 / S2
225               S = SQRT( ONE+TMP*TMP )
226               SESTPR = S2*S
227               C = ( GAMMA / S2 ) / S
228               S = SIGN( ONE, ALPHA ) / S
229            ELSE
230               TMP = S2 / S1
231               C = SQRT( ONE+TMP*TMP )
232               SESTPR = S1*C
233               S = ( ALPHA / S1 ) / C
234               C = SIGN( ONE, GAMMA ) / C
235            END IF
236            RETURN
237         ELSE
238*
239*           normal case
240*
241            ZETA1 = ALPHA / ABSEST
242            ZETA2 = GAMMA / ABSEST
243*
244            B = ( ONE-ZETA1*ZETA1-ZETA2*ZETA2 )*HALF
245            C = ZETA1*ZETA1
246            IF( B.GT.ZERO ) THEN
247               T = C / ( B+SQRT( B*B+C ) )
248            ELSE
249               T = SQRT( B*B+C ) - B
250            END IF
251*
252            SINE = -ZETA1 / T
253            COSINE = -ZETA2 / ( ONE+T )
254            TMP = SQRT( SINE*SINE+COSINE*COSINE )
255            S = SINE / TMP
256            C = COSINE / TMP
257            SESTPR = SQRT( T+ONE )*ABSEST
258            RETURN
259         END IF
260*
261      ELSE IF( JOB.EQ.2 ) THEN
262*
263*        Estimating smallest singular value
264*
265*        special cases
266*
267         IF( SEST.EQ.ZERO ) THEN
268            SESTPR = ZERO
269            IF( MAX( ABSGAM, ABSALP ).EQ.ZERO ) THEN
270               SINE = ONE
271               COSINE = ZERO
272            ELSE
273               SINE = -GAMMA
274               COSINE = ALPHA
275            END IF
276            S1 = MAX( ABS( SINE ), ABS( COSINE ) )
277            S = SINE / S1
278            C = COSINE / S1
279            TMP = SQRT( S*S+C*C )
280            S = S / TMP
281            C = C / TMP
282            RETURN
283         ELSE IF( ABSGAM.LE.EPS*ABSEST ) THEN
284            S = ZERO
285            C = ONE
286            SESTPR = ABSGAM
287            RETURN
288         ELSE IF( ABSALP.LE.EPS*ABSEST ) THEN
289            S1 = ABSGAM
290            S2 = ABSEST
291            IF( S1.LE.S2 ) THEN
292               S = ZERO
293               C = ONE
294               SESTPR = S1
295            ELSE
296               S = ONE
297               C = ZERO
298               SESTPR = S2
299            END IF
300            RETURN
301         ELSE IF( ABSEST.LE.EPS*ABSALP .OR. ABSEST.LE.EPS*ABSGAM ) THEN
302            S1 = ABSGAM
303            S2 = ABSALP
304            IF( S1.LE.S2 ) THEN
305               TMP = S1 / S2
306               C = SQRT( ONE+TMP*TMP )
307               SESTPR = ABSEST*( TMP / C )
308               S = -( GAMMA / S2 ) / C
309               C = SIGN( ONE, ALPHA ) / C
310            ELSE
311               TMP = S2 / S1
312               S = SQRT( ONE+TMP*TMP )
313               SESTPR = ABSEST / S
314               C = ( ALPHA / S1 ) / S
315               S = -SIGN( ONE, GAMMA ) / S
316            END IF
317            RETURN
318         ELSE
319*
320*           normal case
321*
322            ZETA1 = ALPHA / ABSEST
323            ZETA2 = GAMMA / ABSEST
324*
325            NORMA = MAX( ONE+ZETA1*ZETA1+ABS( ZETA1*ZETA2 ),
326     $              ABS( ZETA1*ZETA2 )+ZETA2*ZETA2 )
327*
328*           See if root is closer to zero or to ONE
329*
330            TEST = ONE + TWO*( ZETA1-ZETA2 )*( ZETA1+ZETA2 )
331            IF( TEST.GE.ZERO ) THEN
332*
333*              root is close to zero, compute directly
334*
335               B = ( ZETA1*ZETA1+ZETA2*ZETA2+ONE )*HALF
336               C = ZETA2*ZETA2
337               T = C / ( B+SQRT( ABS( B*B-C ) ) )
338               SINE = ZETA1 / ( ONE-T )
339               COSINE = -ZETA2 / T
340               SESTPR = SQRT( T+FOUR*EPS*EPS*NORMA )*ABSEST
341            ELSE
342*
343*              root is closer to ONE, shift by that amount
344*
345               B = ( ZETA2*ZETA2+ZETA1*ZETA1-ONE )*HALF
346               C = ZETA1*ZETA1
347               IF( B.GE.ZERO ) THEN
348                  T = -C / ( B+SQRT( B*B+C ) )
349               ELSE
350                  T = B - SQRT( B*B+C )
351               END IF
352               SINE = -ZETA1 / T
353               COSINE = -ZETA2 / ( ONE+T )
354               SESTPR = SQRT( ONE+T+FOUR*EPS*EPS*NORMA )*ABSEST
355            END IF
356            TMP = SQRT( SINE*SINE+COSINE*COSINE )
357            S = SINE / TMP
358            C = COSINE / TMP
359            RETURN
360*
361         END IF
362      END IF
363      RETURN
364*
365*     End of SLAIC1
366*
367      END
Note: See TracBrowser for help on using the repository browser.