source: trunk/SOURCES/BLAS/slarft.f @ 23

Last change on this file since 23 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: 10.0 KB
Line 
1*> \brief \b SLARFT
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at 
6*            http://www.netlib.org/lapack/explore-html/ 
7*
8*> \htmlonly
9*> Download SLARFT + dependencies 
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarft.f"> 
11*> [TGZ]</a> 
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarft.f"> 
13*> [ZIP]</a> 
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarft.f"> 
15*> [TXT]</a>
16*> \endhtmlonly 
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
22* 
23*       .. Scalar Arguments ..
24*       CHARACTER          DIRECT, STOREV
25*       INTEGER            K, LDT, LDV, N
26*       ..
27*       .. Array Arguments ..
28*       REAL               T( LDT, * ), TAU( * ), V( LDV, * )
29*       ..
30* 
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> SLARFT forms the triangular factor T of a real block reflector H
38*> of order n, which is defined as a product of k elementary reflectors.
39*>
40*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
41*>
42*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
43*>
44*> If STOREV = 'C', the vector which defines the elementary reflector
45*> H(i) is stored in the i-th column of the array V, and
46*>
47*>    H  =  I - V * T * V**T
48*>
49*> If STOREV = 'R', the vector which defines the elementary reflector
50*> H(i) is stored in the i-th row of the array V, and
51*>
52*>    H  =  I - V**T * T * V
53*> \endverbatim
54*
55*  Arguments:
56*  ==========
57*
58*> \param[in] DIRECT
59*> \verbatim
60*>          DIRECT is CHARACTER*1
61*>          Specifies the order in which the elementary reflectors are
62*>          multiplied to form the block reflector:
63*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
64*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
65*> \endverbatim
66*>
67*> \param[in] STOREV
68*> \verbatim
69*>          STOREV is CHARACTER*1
70*>          Specifies how the vectors which define the elementary
71*>          reflectors are stored (see also Further Details):
72*>          = 'C': columnwise
73*>          = 'R': rowwise
74*> \endverbatim
75*>
76*> \param[in] N
77*> \verbatim
78*>          N is INTEGER
79*>          The order of the block reflector H. N >= 0.
80*> \endverbatim
81*>
82*> \param[in] K
83*> \verbatim
84*>          K is INTEGER
85*>          The order of the triangular factor T (= the number of
86*>          elementary reflectors). K >= 1.
87*> \endverbatim
88*>
89*> \param[in,out] V
90*> \verbatim
91*>          V is REAL array, dimension
92*>                               (LDV,K) if STOREV = 'C'
93*>                               (LDV,N) if STOREV = 'R'
94*>          The matrix V. See further details.
95*> \endverbatim
96*>
97*> \param[in] LDV
98*> \verbatim
99*>          LDV is INTEGER
100*>          The leading dimension of the array V.
101*>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
102*> \endverbatim
103*>
104*> \param[in] TAU
105*> \verbatim
106*>          TAU is REAL array, dimension (K)
107*>          TAU(i) must contain the scalar factor of the elementary
108*>          reflector H(i).
109*> \endverbatim
110*>
111*> \param[out] T
112*> \verbatim
113*>          T is REAL array, dimension (LDT,K)
114*>          The k by k triangular factor T of the block reflector.
115*>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
116*>          lower triangular. The rest of the array is not used.
117*> \endverbatim
118*>
119*> \param[in] LDT
120*> \verbatim
121*>          LDT is INTEGER
122*>          The leading dimension of the array T. LDT >= K.
123*> \endverbatim
124*
125*  Authors:
126*  ========
127*
128*> \author Univ. of Tennessee 
129*> \author Univ. of California Berkeley 
130*> \author Univ. of Colorado Denver 
131*> \author NAG Ltd. 
132*
133*> \date November 2011
134*
135*> \ingroup realOTHERauxiliary
136*
137*> \par Further Details:
138*  =====================
139*>
140*> \verbatim
141*>
142*>  The shape of the matrix V and the storage of the vectors which define
143*>  the H(i) is best illustrated by the following example with n = 5 and
144*>  k = 3. The elements equal to 1 are not stored; the corresponding
145*>  array elements are modified but restored on exit. The rest of the
146*>  array is not used.
147*>
148*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
149*>
150*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
151*>                   ( v1  1    )                     (     1 v2 v2 v2 )
152*>                   ( v1 v2  1 )                     (        1 v3 v3 )
153*>                   ( v1 v2 v3 )
154*>                   ( v1 v2 v3 )
155*>
156*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
157*>
158*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
159*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
160*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
161*>                   (     1 v3 )
162*>                   (        1 )
163*> \endverbatim
164*>
165*  =====================================================================
166      SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
167*
168*  -- LAPACK auxiliary routine (version 3.4.0) --
169*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
170*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171*     November 2011
172*
173*     .. Scalar Arguments ..
174      CHARACTER          DIRECT, STOREV
175      INTEGER            K, LDT, LDV, N
176*     ..
177*     .. Array Arguments ..
178      REAL               T( LDT, * ), TAU( * ), V( LDV, * )
179*     ..
180*
181*  =====================================================================
182*
183*     .. Parameters ..
184      REAL               ONE, ZERO
185      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
186*     ..
187*     .. Local Scalars ..
188      INTEGER            I, J, PREVLASTV, LASTV
189      REAL               VII
190*     ..
191*     .. External Subroutines ..
192      EXTERNAL           SGEMV, STRMV
193*     ..
194*     .. External Functions ..
195      LOGICAL            LSAME
196      EXTERNAL           LSAME
197*     ..
198*     .. Executable Statements ..
199*
200*     Quick return if possible
201*
202      IF( N.EQ.0 )
203     $   RETURN
204*
205      IF( LSAME( DIRECT, 'F' ) ) THEN
206         PREVLASTV = N
207         DO 20 I = 1, K
208            PREVLASTV = MAX( I, PREVLASTV )
209            IF( TAU( I ).EQ.ZERO ) THEN
210*
211*              H(i)  =  I
212*
213               DO 10 J = 1, I
214                  T( J, I ) = ZERO
215   10          CONTINUE
216            ELSE
217*
218*              general case
219*
220               VII = V( I, I )
221               V( I, I ) = ONE
222               IF( LSAME( STOREV, 'C' ) ) THEN
223!                 Skip any trailing zeros.
224                  DO LASTV = N, I+1, -1
225                     IF( V( LASTV, I ).NE.ZERO ) EXIT
226                  END DO
227                  J = MIN( LASTV, PREVLASTV )
228*
229*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
230*
231                  CALL SGEMV( 'Transpose', J-I+1, I-1, -TAU( I ),
232     $                        V( I, 1 ), LDV, V( I, I ), 1, ZERO,
233     $                        T( 1, I ), 1 )
234               ELSE
235!                 Skip any trailing zeros.
236                  DO LASTV = N, I+1, -1
237                     IF( V( I, LASTV ).NE.ZERO ) EXIT
238                  END DO
239                  J = MIN( LASTV, PREVLASTV )
240*
241*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
242*
243                  CALL SGEMV( 'No transpose', I-1, J-I+1, -TAU( I ),
244     $                        V( 1, I ), LDV, V( I, I ), LDV, ZERO,
245     $                        T( 1, I ), 1 )
246               END IF
247               V( I, I ) = VII
248*
249*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
250*
251               CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
252     $                     LDT, T( 1, I ), 1 )
253               T( I, I ) = TAU( I )
254               IF( I.GT.1 ) THEN
255                  PREVLASTV = MAX( PREVLASTV, LASTV )
256               ELSE
257                  PREVLASTV = LASTV
258               END IF
259            END IF
260   20    CONTINUE
261      ELSE
262         PREVLASTV = 1
263         DO 40 I = K, 1, -1
264            IF( TAU( I ).EQ.ZERO ) THEN
265*
266*              H(i)  =  I
267*
268               DO 30 J = I, K
269                  T( J, I ) = ZERO
270   30          CONTINUE
271            ELSE
272*
273*              general case
274*
275               IF( I.LT.K ) THEN
276                  IF( LSAME( STOREV, 'C' ) ) THEN
277                     VII = V( N-K+I, I )
278                     V( N-K+I, I ) = ONE
279!                    Skip any leading zeros.
280                     DO LASTV = 1, I-1
281                        IF( V( LASTV, I ).NE.ZERO ) EXIT
282                     END DO
283                     J = MAX( LASTV, PREVLASTV )
284*
285*                    T(i+1:k,i) :=
286*                            - tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i)
287*
288                     CALL SGEMV( 'Transpose', N-K+I-J+1, K-I, -TAU( I ),
289     $                           V( J, I+1 ), LDV, V( J, I ), 1, ZERO,
290     $                           T( I+1, I ), 1 )
291                     V( N-K+I, I ) = VII
292                  ELSE
293                     VII = V( I, N-K+I )
294                     V( I, N-K+I ) = ONE
295!                    Skip any leading zeros.
296                     DO LASTV = 1, I-1
297                        IF( V( I, LASTV ).NE.ZERO ) EXIT
298                     END DO
299                     J = MAX( LASTV, PREVLASTV )
300*
301*                    T(i+1:k,i) :=
302*                            - tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
303*
304                     CALL SGEMV( 'No transpose', K-I, N-K+I-J+1,
305     $                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
306     $                    ZERO, T( I+1, I ), 1 )
307                     V( I, N-K+I ) = VII
308                  END IF
309*
310*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
311*
312                  CALL STRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
313     $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
314                  IF( I.GT.1 ) THEN
315                     PREVLASTV = MIN( PREVLASTV, LASTV )
316                  ELSE
317                     PREVLASTV = LASTV
318                  END IF
319               END IF
320               T( I, I ) = TAU( I )
321            END IF
322   40    CONTINUE
323      END IF
324      RETURN
325*
326*     End of SLARFT
327*
328      END
Note: See TracBrowser for help on using the repository browser.