1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
SUBROUTINE CLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
*
* -- LAPACK auxiliary routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* September 30, 1994
*
* .. Scalar Arguments ..
CHARACTER SIDE
INTEGER INCV, LDC, M, N
COMPLEX TAU
* ..
* .. Array Arguments ..
COMPLEX C( LDC, * ), V( * ), WORK( * )
* ..
*
* Purpose
* =======
*
* CLARF applies a complex elementary reflector H to a complex M-by-N
* matrix C, from either the left or the right. H is represented in the
* form
*
* H = I - tau * v * v'
*
* where tau is a complex scalar and v is a complex vector.
*
* If tau = 0, then H is taken to be the unit matrix.
*
* To apply H' (the conjugate transpose of H), supply conjg(tau) instead
* tau.
*
* Arguments
* =========
*
* SIDE (input) CHARACTER*1
* = 'L': form H * C
* = 'R': form C * H
*
* M (input) INTEGER
* The number of rows of the matrix C.
*
* N (input) INTEGER
* The number of columns of the matrix C.
*
* V (input) COMPLEX array, dimension
* (1 + (M-1)*abs(INCV)) if SIDE = 'L'
* or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
* The vector v in the representation of H. V is not used if
* TAU = 0.
*
* INCV (input) INTEGER
* The increment between elements of v. INCV <> 0.
*
* TAU (input) COMPLEX
* The value tau in the representation of H.
*
* C (input/output) COMPLEX array, dimension (LDC,N)
* On entry, the M-by-N matrix C.
* On exit, C is overwritten by the matrix H * C if SIDE = 'L',
* or C * H if SIDE = 'R'.
*
* LDC (input) INTEGER
* The leading dimension of the array C. LDC >= max(1,M).
*
* WORK (workspace) COMPLEX array, dimension
* (N) if SIDE = 'L'
* or (M) if SIDE = 'R'
*
* =====================================================================
*
* .. Parameters ..
COMPLEX ONE, ZERO
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ),
$ ZERO = ( 0.0E+0, 0.0E+0 ) )
* ..
* .. External Subroutines ..
EXTERNAL CGEMV, CGERC
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. Executable Statements ..
*
IF( LSAME( SIDE, 'L' ) ) THEN
*
* Form H * C
*
IF( TAU.NE.ZERO ) THEN
*
* w := C' * v
*
CALL CGEMV( 'Conjugate transpose', M, N, ONE, C, LDC, V,
$ INCV, ZERO, WORK, 1 )
*
* C := C - v * w'
*
CALL CGERC( M, N, -TAU, V, INCV, WORK, 1, C, LDC )
END IF
ELSE
*
* Form C * H
*
IF( TAU.NE.ZERO ) THEN
*
* w := C * v
*
CALL CGEMV( 'No transpose', M, N, ONE, C, LDC, V, INCV,
$ ZERO, WORK, 1 )
*
* C := C - w * v'
*
CALL CGERC( M, N, -TAU, WORK, 1, V, INCV, C, LDC )
END IF
END IF
RETURN
*
* End of CLARF
*
END
|
