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
122
123
124
125
|
SUBROUTINE CGERQ2( M, N, A, LDA, TAU, WORK, INFO )
*
* -- LAPACK routine (version 2.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* September 30, 1994
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
* ..
* .. Array Arguments ..
COMPLEX A( LDA, * ), TAU( * ), WORK( * )
* ..
*
* Purpose
* =======
*
* CGERQ2 computes an RQ factorization of a complex m by n matrix A:
* A = R * Q.
*
* Arguments
* =========
*
* M (input) INTEGER
* The number of rows of the matrix A. M >= 0.
*
* N (input) INTEGER
* The number of columns of the matrix A. N >= 0.
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the m by n matrix A.
* On exit, if m <= n, the upper triangle of the subarray
* A(1:m,n-m+1:n) contains the m by m upper triangular matrix R;
* if m >= n, the elements on and above the (m-n)-th subdiagonal
* contain the m by n upper trapezoidal matrix R; the remaining
* elements, with the array TAU, represent the unitary matrix
* Q as a product of elementary reflectors (see Further
* Details).
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
*
* TAU (output) COMPLEX array, dimension (min(M,N))
* The scalar factors of the elementary reflectors (see Further
* Details).
*
* WORK (workspace) COMPLEX array, dimension (M)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* Further Details
* ===============
*
* The matrix Q is represented as a product of elementary reflectors
*
* Q = H(1)' H(2)' . . . H(k)', where k = min(m,n).
*
* Each H(i) has the form
*
* H(i) = I - tau * v * v'
*
* where tau is a complex scalar, and v is a complex vector with
* v(n-k+i+1:n) = 0 and v(n-k+i) = 1; conjg(v(1:n-k+i-1)) is stored on
* exit in A(m-k+i,1:n-k+i-1), and tau in TAU(i).
*
* =====================================================================
*
* .. Parameters ..
COMPLEX ONE
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
INTEGER I, K
COMPLEX ALPHA
* ..
* .. External Subroutines ..
EXTERNAL CLACGV, CLARF, CLARFG, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input arguments
*
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'CGERQ2', -INFO )
RETURN
END IF
*
K = MIN( M, N )
*
DO 10 I = K, 1, -1
*
* Generate elementary reflector H(i) to annihilate
* A(m-k+i,1:n-k+i-1)
*
CALL CLACGV( N-K+I, A( M-K+I, 1 ), LDA )
ALPHA = A( M-K+I, N-K+I )
CALL CLARFG( N-K+I, ALPHA, A( M-K+I, 1 ), LDA,
$ TAU( I ) )
*
* Apply H(i) to A(1:m-k+i-1,1:n-k+i) from the right
*
A( M-K+I, N-K+I ) = ONE
CALL CLARF( 'Right', M-K+I-1, N-K+I, A( M-K+I, 1 ), LDA,
$ TAU( I ), A, LDA, WORK )
A( M-K+I, N-K+I ) = ALPHA
CALL CLACGV( N-K+I-1, A( M-K+I, 1 ), LDA )
10 CONTINUE
RETURN
*
* End of CGERQ2
*
END
|
