Skip to content

Commit 6c0a98f

Browse files
implement clarf1f, #1011
1 parent 8dd7e13 commit 6c0a98f

File tree

3 files changed

+269
-2
lines changed

3 files changed

+269
-2
lines changed

SRC/CMakeLists.txt

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -218,7 +218,7 @@ set(CLASRC
218218
claqhb.f claqhe.f claqhp.f claqp2.f claqps.f claqp2rk.f claqp3rk.f claqsb.f
219219
claqr0.f claqr1.f claqr2.f claqr3.f claqr4.f claqr5.f
220220
claqsp.f claqsy.f clar1v.f clar2v.f ilaclr.f ilaclc.f
221-
clarf.f clarfb.f clarfb_gett.f clarfg.f clarfgp.f clarft.f
221+
clarf.f clarf1f.f clarfb.f clarfb_gett.f clarfg.f clarfgp.f clarft.f
222222
clarfx.f clarfy.f clargv.f clarnv.f clarrv.f clartg.f90 clartv.f
223223
clarz.f clarzb.f clarzt.f clascl.f claset.f clasr.f classq.f90
224224
claswp.f clasyf.f clasyf_rook.f clasyf_rk.f clasyf_aa.f

SRC/Makefile

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -249,7 +249,7 @@ CLASRC = \
249249
claqhb.o claqhe.o claqhp.o claqp2.o claqps.o claqp2rk.o claqp3rk.o claqsb.o \
250250
claqr0.o claqr1.o claqr2.o claqr3.o claqr4.o claqr5.o \
251251
claqsp.o claqsy.o clar1v.o clar2v.o ilaclr.o ilaclc.o \
252-
clarf.o clarfb.o clarfb_gett.o clarfg.o clarft.o clarfgp.o \
252+
clarf.o clarf1f.o clarfb.o clarfb_gett.o clarfg.o clarft.o clarfgp.o \
253253
clarfx.o clarfy.o clargv.o clarnv.o clarrv.o clartg.o clartv.o \
254254
clarz.o clarzb.o clarzt.o clascl.o claset.o clasr.o classq.o \
255255
claswp.o clasyf.o clasyf_rook.o clasyf_rk.o clasyf_aa.o \

SRC/clarf1f.f

Lines changed: 267 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,267 @@
1+
*> \brief \b CLARF1F applies an elementary reflector to a general rectangular
2+
* matrix assuming v(1) = 1.
3+
*
4+
* =========== DOCUMENTATION ===========
5+
*
6+
* Online html documentation available at
7+
* http://www.netlib.org/lapack/explore-html/
8+
*
9+
*> \htmlonly
10+
*> Download CLARF1F + dependencies
11+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarf.f">
12+
*> [TGZ]</a>
13+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarf.f">
14+
*> [ZIP]</a>
15+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarf.f">
16+
*> [TXT]</a>
17+
*> \endhtmlonly
18+
*
19+
* Definition:
20+
* ===========
21+
*
22+
* SUBROUTINE CLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
23+
*
24+
* .. Scalar Arguments ..
25+
* CHARACTER SIDE
26+
* INTEGER INCV, LDC, M, N
27+
* COMPLEX TAU
28+
* ..
29+
* .. Array Arguments ..
30+
* COMPLEX C( LDC, * ), V( * ), WORK( * )
31+
* ..
32+
*
33+
*
34+
*> \par Purpose:
35+
* =============
36+
*>
37+
*> \verbatim
38+
*>
39+
*> CLARF1F applies a complex elementary reflector H to a complex m by n matrix
40+
*> C, from either the left or the right. H is represented in the form
41+
*>
42+
*> H = I - tau * v * v**H
43+
*>
44+
*> where tau is a complex scalar and v is a complex vector assuming v(1) = 1.
45+
*>
46+
*> If tau = 0, then H is taken to be the unit matrix.
47+
*>
48+
*> To apply H**H (the conjugate transpose of H), supply conjg(tau) instead
49+
*> tau.
50+
*> \endverbatim
51+
*
52+
* Arguments:
53+
* ==========
54+
*
55+
*> \param[in] SIDE
56+
*> \verbatim
57+
*> SIDE is CHARACTER*1
58+
*> = 'L': form H * C
59+
*> = 'R': form C * H
60+
*> \endverbatim
61+
*>
62+
*> \param[in] M
63+
*> \verbatim
64+
*> M is INTEGER
65+
*> The number of rows of the matrix C.
66+
*> \endverbatim
67+
*>
68+
*> \param[in] N
69+
*> \verbatim
70+
*> N is INTEGER
71+
*> The number of columns of the matrix C.
72+
*> \endverbatim
73+
*>
74+
*> \param[in] V
75+
*> \verbatim
76+
*> V is COMPLEX array, dimension
77+
*> (1 + (M-1)*abs(INCV)) if SIDE = 'L'
78+
*> or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
79+
*> The vector v in the representation of H. V is not used if
80+
*> TAU = 0.
81+
*> \endverbatim
82+
*>
83+
*> \param[in] INCV
84+
*> \verbatim
85+
*> INCV is INTEGER
86+
*> The increment between elements of v. INCV <> 0.
87+
*> \endverbatim
88+
*>
89+
*> \param[in] TAU
90+
*> \verbatim
91+
*> TAU is COMPLEX
92+
*> The value tau in the representation of H.
93+
*> \endverbatim
94+
*>
95+
*> \param[in,out] C
96+
*> \verbatim
97+
*> C is COMPLEX array, dimension (LDC,N)
98+
*> On entry, the m by n matrix C.
99+
*> On exit, C is overwritten by the matrix H * C if SIDE = 'L',
100+
*> or C * H if SIDE = 'R'.
101+
*> \endverbatim
102+
*>
103+
*> \param[in] LDC
104+
*> \verbatim
105+
*> LDC is INTEGER
106+
*> The leading dimension of the array C. LDC >= max(1,M).
107+
*> \endverbatim
108+
*>
109+
*> \param[out] WORK
110+
*> \verbatim
111+
*> WORK is COMPLEX array, dimension
112+
*> (N) if SIDE = 'L'
113+
*> or (M) if SIDE = 'R'
114+
*> \endverbatim
115+
*
116+
* Authors:
117+
* ========
118+
*
119+
*> \author Univ. of Tennessee
120+
*> \author Univ. of California Berkeley
121+
*> \author Univ. of Colorado Denver
122+
*> \author NAG Ltd.
123+
*
124+
*> \ingroup larf1f
125+
*
126+
* =====================================================================
127+
SUBROUTINE CLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
128+
*
129+
* -- LAPACK auxiliary routine --
130+
* -- LAPACK is a software package provided by Univ. of Tennessee, --
131+
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
132+
*
133+
* .. Scalar Arguments ..
134+
CHARACTER SIDE
135+
INTEGER INCV, LDC, M, N
136+
COMPLEX TAU
137+
* ..
138+
* .. Array Arguments ..
139+
COMPLEX C( LDC, * ), V( * ), WORK( * )
140+
* ..
141+
*
142+
* =====================================================================
143+
*
144+
* .. Parameters ..
145+
COMPLEX ONE, ZERO
146+
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ),
147+
$ ZERO = ( 0.0E+0, 0.0E+0 ) )
148+
* ..
149+
* .. Local Scalars ..
150+
LOGICAL APPLYLEFT
151+
INTEGER I, LASTV, LASTC
152+
* ..
153+
* .. External Subroutines ..
154+
EXTERNAL CGEMV, CGER, CSCAL
155+
* ..
156+
* .. Intrinsic Functions ..
157+
INTRINSIC CONJG
158+
* ..
159+
* .. External Functions ..
160+
LOGICAL LSAME
161+
INTEGER ILACLR, ILACLC
162+
EXTERNAL LSAME, ILACLR, ILACLC
163+
* ..
164+
* .. Executable Statements ..
165+
*
166+
APPLYLEFT = LSAME( SIDE, 'L' )
167+
LASTV = 0
168+
LASTC = 0
169+
IF( TAU.NE.ZERO ) THEN
170+
! Set up variables for scanning V. LASTV begins pointing to the end
171+
! of V.
172+
IF( APPLYLEFT ) THEN
173+
LASTV = M
174+
ELSE
175+
LASTV = N
176+
END IF
177+
IF( INCV.GT.0 ) THEN
178+
I = 1 + (LASTV-1) * INCV
179+
ELSE
180+
I = 1
181+
END IF
182+
! Look for the last non-zero row in V.
183+
DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
184+
LASTV = LASTV - 1
185+
I = I - INCV
186+
END DO
187+
IF( APPLYLEFT ) THEN
188+
! Scan for the last non-zero column in C(1:lastv,:).
189+
LASTC = ILACLC(LASTV, N, C, LDC)
190+
ELSE
191+
! Scan for the last non-zero row in C(:,1:lastv).
192+
LASTC = ILACLR(M, LASTV, C, LDC)
193+
END IF
194+
END IF
195+
IF( LASTC.EQ.0 .OR. LASTV.EQ.0 ) THEN
196+
RETURN
197+
END IF
198+
IF( APPLYLEFT ) THEN
199+
*
200+
* Form H * C
201+
*
202+
IF( LASTV.EQ.1 ) THEN
203+
*
204+
* C(1,1:lastc) := ( 1 - tau ) * C(1,1:lastc)
205+
*
206+
CALL CSCAL( LASTC, ONE - TAU, C, LDC )
207+
ELSE
208+
*
209+
* w(1:lastc,1) := C(2:lastv,1:lastc)**H * v(2:lastv,1)
210+
*
211+
CALL CGEMV( 'Conjugate transpose', LASTV - 1, LASTC, ONE,
212+
$ C( 2, 1 ), LDC, V( 1 + INCV ), INCV, ZERO,
213+
$ WORK, 1 )
214+
*
215+
* w(1:lastc,1) += v(1,1) * C(1,1:lastc)**H
216+
*
217+
DO I = 1, LASTC
218+
WORK( I ) = WORK( I ) + CONJG( C( 1, I ) )
219+
END DO
220+
*
221+
* C(1, 1:lastc) += - tau * v(1,1) * w(1:lastc,1)**H
222+
*
223+
DO I = 1, LASTC
224+
C( 1, I ) = C( 1, I ) - TAU * CONJG( WORK( I ) )
225+
END DO
226+
*
227+
* C(2:lastv,1:lastc) += - tau * v(2:lastv,1) * w(1:lastc,1)**H
228+
*
229+
CALL CGERC( LASTV - 1, LASTC, -TAU, V( 1 + INCV ), INCV,
230+
$ WORK, 1, C( 2, 1 ), LDC )
231+
END IF
232+
ELSE
233+
*
234+
* Form C * H
235+
*
236+
IF( LASTV.EQ.1 ) THEN
237+
*
238+
* C(1:lastc,1) := ( 1 - tau ) * C(1:lastc,1)
239+
*
240+
CALL CSCAL( LASTC, ONE - TAU, C, 1 )
241+
ELSE
242+
*
243+
* w(1:lastc,1) := C(1:lastc,2:lastv) * v(2:lastv,1)
244+
*
245+
CALL CGEMV( 'No transpose', LASTC, LASTV - 1, ONE,
246+
$ C( 1, 2 ), LDC, V( 1 + INCV ), INCV, ZERO,
247+
$ WORK, 1 )
248+
*
249+
* w(1:lastc,1) += v(1,1) * C(1:lastc,1)
250+
*
251+
CALL CAXPY( LASTC, ONE, C, 1, WORK, 1 )
252+
*
253+
* C(1:lastc,1) += - tau * v(1,1) * w(1:lastc,1)
254+
*
255+
CALL CAXPY( LASTC, -TAU, WORK, 1, C, 1 )
256+
*
257+
* C(1:lastc,2:lastv) += - tau * w(1:lastc,1) * v(2:lastv)**H
258+
*
259+
CALL CGERC( LASTC, LASTV - 1, -TAU, WORK, 1,
260+
$ V( 1 + INCV ), INCV, C( 1, 2 ), LDC )
261+
END IF
262+
END IF
263+
RETURN
264+
*
265+
* End of CLARF1F
266+
*
267+
END

0 commit comments

Comments
 (0)