Skip to content

Commit 7708f1e

Browse files
jprhynejprhyne
committed
update dlarf1f.f and zlarf1f.f to not reference v(1)
1 parent 468cb59 commit 7708f1e

File tree

3 files changed

+116
-54
lines changed

3 files changed

+116
-54
lines changed

‎SRC/dlarf1f.f

Lines changed: 77 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,7 @@
7474
*> (1 + (M-1)*abs(INCV)) if SIDE = 'L'
7575
*> or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
7676
*> The vector v in the representation of H. V is not used if
77-
*> TAU = 0.
77+
*> TAU = 0. V(1) is not referenced or modified.
7878
*> \endverbatim
7979
*>
8080
*> \param[in] INCV
@@ -110,6 +110,40 @@
110110
*> or (M) if SIDE = 'R'
111111
*> \endverbatim
112112
*
113+
* To take advantage of the fact that v(1) = 1, we do the following
114+
* v = [ 1 v_2 ]**T
115+
* If SIDE='L'
116+
* |-----|
117+
* | C_1 |
118+
* C =| C_2 |
119+
* |-----|
120+
* C_1\in\mathbb{R}^{1\times n}, C_2\in\mathbb{R}^{m-1\times n}
121+
* So we compute:
122+
* C = HC = (I - \tau vv**T)C
123+
* = C - \tau vv**T C
124+
* w = C**T v = [ C_1**T C_2**T ] [ 1 v_2 ]**T
125+
* = C_1**T + C_2**T v ( DGEMM then DAXPY )
126+
* C = C - \tau vv**T C
127+
* = C - \tau vw**T
128+
* Giving us C_1 = C_1 - \tau w**T ( DAXPY )
129+
* and
130+
* C_2 = C_2 - \tau v_2w**T ( DGER )
131+
* If SIDE='R'
132+
*
133+
* C = [ C_1 C_2 ]
134+
* C_1\in\mathbb{R}^{m\times 1}, C_2\in\mathbb{R}^{m\times n-1}
135+
* So we compute:
136+
* C = CH = C(I - \tau vv**T)
137+
* = C - \tau Cvv**T
138+
*
139+
* w = Cv = [ C_1 C_2 ] [ 1 v_2 ]**T
140+
* = C_1 + C_2v_2 ( DGEMM then DAXPY )
141+
* C = C - \tau Cvv**T
142+
* = C - \tau wv**T
143+
* Giving us C_1 = C_1 - \tau w ( DAXPY )
144+
* and
145+
* C_2 = C_2 - \tau wv_2**T ( DGER )
146+
*
113147
* Authors:
114148
* ========
115149
*
@@ -175,7 +209,9 @@ SUBROUTINE DLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
175209
I = 1
176210
END IF
177211
! Look for the last non-zero row in V.
178-
DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
212+
! Since we are assuming that V(1) = 1, and it is not stored, so we
213+
! shouldn't access it.
214+
DO WHILE( LASTV.GT.1 .AND. V( I ).EQ.ZERO )
179215
LASTV = LASTV - 1
180216
I = I - INCV
181217
END DO
@@ -186,67 +222,63 @@ SUBROUTINE DLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
186222
! Scan for the last non-zero row in C(:,1:lastv).
187223
LASTC = ILADLR(M, LASTV, C, LDC)
188224
END IF
189-
END IF
190-
IF( LASTC.EQ.0 .OR. LASTV.EQ.0 ) THEN
225+
ELSE
226+
! TAU is 0, so H = I. Meaning HC = C = CH.
191227
RETURN
192228
END IF
193229
IF( APPLYLEFT ) THEN
194230
*
195231
* Form H * C
196232
*
197-
IF( LASTV.GT.0 ) THEN
198-
! Check if m = 1. This means v = 1, So we just need to compute
199-
! C := HC = (1-\tau)C.
200-
IF( M.EQ.1 .OR. LASTV.EQ.1) THEN
201-
CALL DSCAL(LASTC, ONE - TAU, C, LDC)
202-
ELSE
233+
! Check if lastv = 1. This means v = 1, So we just need to compute
234+
! C := HC = (1-\tau)C.
235+
IF( LASTV.EQ.1 ) THEN
236+
CALL DSCAL(LASTC, ONE - TAU, C, LDC)
237+
ELSE
203238
*
204-
* w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1)
239+
* w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1)
205240
*
206-
! w(1:lastc,1) := C(2:lastv,1:lastc)**T * v(2:lastv,1)
207-
CALL DGEMV( 'Transpose', LASTV-1, LASTC, ONE, C(1+1,1),
208-
$ LDC, V(1+INCV), INCV, ZERO, WORK, 1)
209-
! w(1:lastc,1) += C(1,1:lastc)**T * v(1,1) = C(1,1:lastc)**T
210-
CALL DAXPY(LASTC, ONE, C, LDC, WORK, 1)
241+
! w(1:lastc,1) := C(2:lastv,1:lastc)**T * v(2:lastv,1)
242+
CALL DGEMV( 'Transpose', LASTV-1, LASTC, ONE, C(1+1,1),
243+
$ LDC, V(1+INCV), INCV, ZERO, WORK, 1)
244+
! w(1:lastc,1) += C(1,1:lastc)**T * v(1,1) = C(1,1:lastc)**T
245+
CALL DAXPY(LASTC, ONE, C, LDC, WORK, 1)
211246
*
212247
* C(1:lastv,1:lastc) := C(...) - tau * v(1:lastv,1) * w(1:lastc,1)**T
213248
*
214249
! C(1, 1:lastc) := C(...) - tau * v(1,1) * w(1:lastc,1)**T
215250
! = C(...) - tau * w(1:lastc,1)**T
216-
CALL DAXPY(LASTC, -TAU, WORK, 1, C, LDC)
217-
! C(2:lastv,1:lastc) := C(...) - tau * v(2:lastv,1)*w(1:lastc,1)**T
218-
CALL DGER(LASTV-1, LASTC, -TAU, V(1+INCV), INCV, WORK, 1,
219-
$ C(1+1,1), LDC)
220-
END IF
251+
CALL DAXPY(LASTC, -TAU, WORK, 1, C, LDC)
252+
! C(2:lastv,1:lastc) := C(...) - tau * v(2:lastv,1)*w(1:lastc,1)**T
253+
CALL DGER(LASTV-1, LASTC, -TAU, V(1+INCV), INCV, WORK, 1,
254+
$ C(1+1,1), LDC)
221255
END IF
222256
ELSE
223257
*
224258
* Form C * H
225259
*
226-
IF( LASTV.GT.0 ) THEN
227-
! Check if n = 1. This means v = 1, so we just need to compute
228-
! C := CH = C(1-\tau).
229-
IF( N.EQ.1 .OR. LASTV.EQ.1) THEN
230-
CALL DSCAL(LASTC, ONE - TAU, C, 1)
231-
ELSE
232-
*
233-
* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)
234-
*
235-
! w(1:lastc,1) := C(1:lastc,2:lastv) * v(2:lastv,1)
236-
CALL DGEMV( 'No transpose', LASTC, LASTV-1, ONE,
237-
$ C(1,1+1), LDC, V(1+INCV), INCV, ZERO, WORK, 1 )
238-
! w(1:lastc,1) += C(1:lastc,1) v(1,1) = C(1:lastc,1)
239-
CALL DAXPY(LASTC, ONE, C, 1, WORK, 1)
240-
*
241-
* C(1:lastc,1:lastv) := C(...) - tau * w(1:lastc,1) * v(1:lastv,1)**T
242-
*
243-
! C(1:lastc,1) := C(...) - tau * w(1:lastc,1) * v(1,1)**T
244-
! = C(...) - tau * w(1:lastc,1)
245-
CALL DAXPY(LASTC, -TAU, WORK, 1, C, 1)
246-
! C(1:lastc,2:lastv) := C(...) - tau * w(1:lastc,1) * v(2:lastv)**T
247-
CALL DGER( LASTC, LASTV-1, -TAU, WORK, 1, V(1+INCV),
248-
$ INCV, C(1,1+1), LDC )
249-
END IF
260+
! Check if n = 1. This means v = 1, so we just need to compute
261+
! C := CH = C(1-\tau).
262+
IF( LASTV.EQ.1 ) THEN
263+
CALL DSCAL(LASTC, ONE - TAU, C, 1)
264+
ELSE
265+
*
266+
* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)
267+
*
268+
! w(1:lastc,1) := C(1:lastc,2:lastv) * v(2:lastv,1)
269+
CALL DGEMV( 'No transpose', LASTC, LASTV-1, ONE,
270+
$ C(1,1+1), LDC, V(1+INCV), INCV, ZERO, WORK, 1 )
271+
! w(1:lastc,1) += C(1:lastc,1) v(1,1) = C(1:lastc,1)
272+
CALL DAXPY(LASTC, ONE, C, 1, WORK, 1)
273+
*
274+
* C(1:lastc,1:lastv) := C(...) - tau * w(1:lastc,1) * v(1:lastv,1)**T
275+
*
276+
! C(1:lastc,1) := C(...) - tau * w(1:lastc,1) * v(1,1)**T
277+
! = C(...) - tau * w(1:lastc,1)
278+
CALL DAXPY(LASTC, -TAU, WORK, 1, C, 1)
279+
! C(1:lastc,2:lastv) := C(...) - tau * w(1:lastc,1) * v(2:lastv)**T
280+
CALL DGER( LASTC, LASTV-1, -TAU, WORK, 1, V(1+INCV),
281+
$ INCV, C(1,1+1), LDC )
250282
END IF
251283
END IF
252284
RETURN

‎SRC/dorg2l.f

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -176,7 +176,6 @@ SUBROUTINE DORG2L( M, N, K, A, LDA, TAU, WORK, INFO )
176176
*
177177
* Apply H(i) to A(1:m-k+i,1:n-k+i) from the left
178178
*
179-
A( M-N+II, II ) = ONE
180179
CALL DLARF1L( 'Left', M-N+II, II-1, A( 1, II ), 1, TAU( I ),
181180
$ A,
182181
$ LDA, WORK )

‎SRC/zlarf1f.f

Lines changed: 39 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,7 @@
7575
*> (1 + (M-1)*abs(INCV)) if SIDE = 'L'
7676
*> or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
7777
*> The vector v in the representation of H. V is not used if
78-
*> TAU = 0.
78+
*> TAU = 0. V(1) is not referenced or modified.
7979
*> \endverbatim
8080
*>
8181
*> \param[in] INCV
@@ -110,6 +110,39 @@
110110
*> (N) if SIDE = 'L'
111111
*> or (M) if SIDE = 'R'
112112
*> \endverbatim
113+
* To take advantage of the fact that v(1) = 1, we do the following
114+
* v = [ 1 v_2 ]**T
115+
* If SIDE='L'
116+
* |-----|
117+
* | C_1 |
118+
* C =| C_2 |
119+
* |-----|
120+
* C_1\in\mathbb{C}^{1\times n}, C_2\in\mathbb{C}^{m-1\times n}
121+
* So we compute:
122+
* C = HC = (I - \tau vv**T)C
123+
* = C - \tau vv**T C
124+
* w = C**T v = [ C_1**T C_2**T ] [ 1 v_2 ]**T
125+
* = C_1**T + C_2**T v ( ZGEMM then ZAXPYC-like )
126+
* C = C - \tau vv**T C
127+
* = C - \tau vw**T
128+
* Giving us C_1 = C_1 - \tau w**T ( ZAXPYC-like )
129+
* and
130+
* C_2 = C_2 - \tau v_2w**T ( ZGERC )
131+
* If SIDE='R'
132+
*
133+
* C = [ C_1 C_2 ]
134+
* C_1\in\mathbb{C}^{m\times 1}, C_2\in\mathbb{C}^{m\times n-1}
135+
* So we compute:
136+
* C = CH = C(I - \tau vv**T)
137+
* = C - \tau Cvv**T
138+
*
139+
* w = Cv = [ C_1 C_2 ] [ 1 v_2 ]**T
140+
* = C_1 + C_2v_2 ( ZGEMM then ZAXPYC-like )
141+
* C = C - \tau Cvv**T
142+
* = C - \tau wv**T
143+
* Giving us C_1 = C_1 - \tau w ( ZAXPYC-like )
144+
* and
145+
* C_2 = C_2 - \tau wv_2**T ( ZGERC )
113146
*
114147
* Authors:
115148
* ========
@@ -177,7 +210,9 @@ SUBROUTINE ZLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
177210
I = 1
178211
END IF
179212
! Look for the last non-zero row in V.
180-
DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
213+
! Since we are assuming that V(1) = 1, and it is not stored, so we
214+
! shouldn't access it.
215+
DO WHILE( LASTV.GT.1 .AND. V( I ).EQ.ZERO )
181216
LASTV = LASTV - 1
182217
I = I - INCV
183218
END DO
@@ -196,10 +231,9 @@ SUBROUTINE ZLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
196231
*
197232
* Form H * C
198233
*
199-
IF( LASTV.GT.0 ) THEN
200234
! Check if m = 1. This means v = 1, So we just need to compute
201235
! C := HC = (1-\tau)C.
202-
IF( M.EQ.1 .OR. LASTV.EQ.1) THEN
236+
IF( LASTV.EQ.1 ) THEN
203237
CALL ZSCAL(LASTC, ONE - TAU, C, LDC)
204238
ELSE
205239
*
@@ -230,15 +264,13 @@ SUBROUTINE ZLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
230264
CALL ZGERC(LASTV-1, LASTC, -TAU, V(1+INCV), INCV, WORK,
231265
$ 1, C(1+1,1), LDC)
232266
END IF
233-
END IF
234267
ELSE
235268
*
236269
* Form C * H
237270
*
238-
IF( LASTV.GT.0 ) THEN
239271
! Check if n = 1. This means v = 1, so we just need to compute
240272
! C := CH = C(1-\tau).
241-
IF( N.EQ.1 .OR. LASTV.EQ.1) THEN
273+
IF( LASTV.EQ.1 ) THEN
242274
CALL ZSCAL(LASTC, ONE - TAU, C, 1)
243275
ELSE
244276
*
@@ -259,7 +291,6 @@ SUBROUTINE ZLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
259291
CALL ZGERC( LASTC, LASTV-1, -TAU, WORK, 1, V(1+INCV),
260292
$ INCV, C(1,1+1), LDC )
261293
END IF
262-
END IF
263294
END IF
264295
RETURN
265296
*

0 commit comments

Comments
 (0)