Skip to content

Commit b8b9771

Browse files
EduardFedorenkovEduardFedorenkov
committed
implement clarf1l, #1011 
1 parent ea943fc commit b8b9771

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 clarf1f.f clarfb.f clarfb_gett.f clarfg.f clarfgp.f clarft.f
221+
clarf.f clarf1f.f clarf1l.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 clarf1f.o clarfb.o clarfb_gett.o clarfg.o clarft.o clarfgp.o \
252+
clarf.o clarf1f.o clarf1l.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/clarf1l.f

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

0 commit comments

Comments
 (0)