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