Annotation of rpl/lapack/lapack/dgerqf.f, revision 1.6
1.1 bertrand 1: SUBROUTINE DGERQF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, LDA, LWORK, M, N
10: * ..
11: * .. Array Arguments ..
12: DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
13: * ..
14: *
15: * Purpose
16: * =======
17: *
18: * DGERQF computes an RQ factorization of a real M-by-N matrix A:
19: * A = R * Q.
20: *
21: * Arguments
22: * =========
23: *
24: * M (input) INTEGER
25: * The number of rows of the matrix A. M >= 0.
26: *
27: * N (input) INTEGER
28: * The number of columns of the matrix A. N >= 0.
29: *
30: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
31: * On entry, the M-by-N matrix A.
32: * On exit,
33: * if m <= n, the upper triangle of the subarray
34: * A(1:m,n-m+1:n) contains the M-by-M upper triangular matrix R;
35: * if m >= n, the elements on and above the (m-n)-th subdiagonal
36: * contain the M-by-N upper trapezoidal matrix R;
37: * the remaining elements, with the array TAU, represent the
38: * orthogonal matrix Q as a product of min(m,n) elementary
39: * reflectors (see Further Details).
40: *
41: * LDA (input) INTEGER
42: * The leading dimension of the array A. LDA >= max(1,M).
43: *
44: * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
45: * The scalar factors of the elementary reflectors (see Further
46: * Details).
47: *
48: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
49: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
50: *
51: * LWORK (input) INTEGER
52: * The dimension of the array WORK. LWORK >= max(1,M).
53: * For optimum performance LWORK >= M*NB, where NB is
54: * the optimal blocksize.
55: *
56: * If LWORK = -1, then a workspace query is assumed; the routine
57: * only calculates the optimal size of the WORK array, returns
58: * this value as the first entry of the WORK array, and no error
59: * message related to LWORK is issued by XERBLA.
60: *
61: * INFO (output) INTEGER
62: * = 0: successful exit
63: * < 0: if INFO = -i, the i-th argument had an illegal value
64: *
65: * Further Details
66: * ===============
67: *
68: * The matrix Q is represented as a product of elementary reflectors
69: *
70: * Q = H(1) H(2) . . . H(k), where k = min(m,n).
71: *
72: * Each H(i) has the form
73: *
74: * H(i) = I - tau * v * v'
75: *
76: * where tau is a real scalar, and v is a real vector with
77: * v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in
78: * A(m-k+i,1:n-k+i-1), and tau in TAU(i).
79: *
80: * =====================================================================
81: *
82: * .. Local Scalars ..
83: LOGICAL LQUERY
84: INTEGER I, IB, IINFO, IWS, K, KI, KK, LDWORK, LWKOPT,
85: $ MU, NB, NBMIN, NU, NX
86: * ..
87: * .. External Subroutines ..
88: EXTERNAL DGERQ2, DLARFB, DLARFT, XERBLA
89: * ..
90: * .. Intrinsic Functions ..
91: INTRINSIC MAX, MIN
92: * ..
93: * .. External Functions ..
94: INTEGER ILAENV
95: EXTERNAL ILAENV
96: * ..
97: * .. Executable Statements ..
98: *
99: * Test the input arguments
100: *
101: INFO = 0
102: LQUERY = ( LWORK.EQ.-1 )
103: IF( M.LT.0 ) THEN
104: INFO = -1
105: ELSE IF( N.LT.0 ) THEN
106: INFO = -2
107: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
108: INFO = -4
109: END IF
110: *
111: IF( INFO.EQ.0 ) THEN
112: K = MIN( M, N )
113: IF( K.EQ.0 ) THEN
114: LWKOPT = 1
115: ELSE
116: NB = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
117: LWKOPT = M*NB
118: END IF
119: WORK( 1 ) = LWKOPT
120: *
121: IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
122: INFO = -7
123: END IF
124: END IF
125: *
126: IF( INFO.NE.0 ) THEN
127: CALL XERBLA( 'DGERQF', -INFO )
128: RETURN
129: ELSE IF( LQUERY ) THEN
130: RETURN
131: END IF
132: *
133: * Quick return if possible
134: *
135: IF( K.EQ.0 ) THEN
136: RETURN
137: END IF
138: *
139: NBMIN = 2
140: NX = 1
141: IWS = M
142: IF( NB.GT.1 .AND. NB.LT.K ) THEN
143: *
144: * Determine when to cross over from blocked to unblocked code.
145: *
146: NX = MAX( 0, ILAENV( 3, 'DGERQF', ' ', M, N, -1, -1 ) )
147: IF( NX.LT.K ) THEN
148: *
149: * Determine if workspace is large enough for blocked code.
150: *
151: LDWORK = M
152: IWS = LDWORK*NB
153: IF( LWORK.LT.IWS ) THEN
154: *
155: * Not enough workspace to use optimal NB: reduce NB and
156: * determine the minimum value of NB.
157: *
158: NB = LWORK / LDWORK
159: NBMIN = MAX( 2, ILAENV( 2, 'DGERQF', ' ', M, N, -1,
160: $ -1 ) )
161: END IF
162: END IF
163: END IF
164: *
165: IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
166: *
167: * Use blocked code initially.
168: * The last kk rows are handled by the block method.
169: *
170: KI = ( ( K-NX-1 ) / NB )*NB
171: KK = MIN( K, KI+NB )
172: *
173: DO 10 I = K - KK + KI + 1, K - KK + 1, -NB
174: IB = MIN( K-I+1, NB )
175: *
176: * Compute the RQ factorization of the current block
177: * A(m-k+i:m-k+i+ib-1,1:n-k+i+ib-1)
178: *
179: CALL DGERQ2( IB, N-K+I+IB-1, A( M-K+I, 1 ), LDA, TAU( I ),
180: $ WORK, IINFO )
181: IF( M-K+I.GT.1 ) THEN
182: *
183: * Form the triangular factor of the block reflector
184: * H = H(i+ib-1) . . . H(i+1) H(i)
185: *
186: CALL DLARFT( 'Backward', 'Rowwise', N-K+I+IB-1, IB,
187: $ A( M-K+I, 1 ), LDA, TAU( I ), WORK, LDWORK )
188: *
189: * Apply H to A(1:m-k+i-1,1:n-k+i+ib-1) from the right
190: *
191: CALL DLARFB( 'Right', 'No transpose', 'Backward',
192: $ 'Rowwise', M-K+I-1, N-K+I+IB-1, IB,
193: $ A( M-K+I, 1 ), LDA, WORK, LDWORK, A, LDA,
194: $ WORK( IB+1 ), LDWORK )
195: END IF
196: 10 CONTINUE
197: MU = M - K + I + NB - 1
198: NU = N - K + I + NB - 1
199: ELSE
200: MU = M
201: NU = N
202: END IF
203: *
204: * Use unblocked code to factor the last or only block
205: *
206: IF( MU.GT.0 .AND. NU.GT.0 )
207: $ CALL DGERQ2( MU, NU, A, LDA, TAU, WORK, IINFO )
208: *
209: WORK( 1 ) = IWS
210: RETURN
211: *
212: * End of DGERQF
213: *
214: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgerqf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack