Annotation of rpl/lapack/lapack/dgeqrf.f, revision 1.18
1.9 bertrand 1: *> \brief \b DGEQRF
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.15 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.15 bertrand 9: *> Download DGEQRF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqrf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqrf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqrf.f">
1.9 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
1.15 bertrand 22: *
1.9 bertrand 23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, LWORK, M, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
28: * ..
1.15 bertrand 29: *
1.9 bertrand 30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DGEQRF computes a QR factorization of a real M-by-N matrix A:
1.18 ! bertrand 37: *>
! 38: *> A = Q * ( R ),
! 39: *> ( 0 )
! 40: *>
! 41: *> where:
! 42: *>
! 43: *> Q is a M-by-M orthogonal matrix;
! 44: *> R is an upper-triangular N-by-N matrix;
! 45: *> 0 is a (M-N)-by-N zero matrix, if M > N.
! 46: *>
1.9 bertrand 47: *> \endverbatim
48: *
49: * Arguments:
50: * ==========
51: *
52: *> \param[in] M
53: *> \verbatim
54: *> M is INTEGER
55: *> The number of rows of the matrix A. M >= 0.
56: *> \endverbatim
57: *>
58: *> \param[in] N
59: *> \verbatim
60: *> N is INTEGER
61: *> The number of columns of the matrix A. N >= 0.
62: *> \endverbatim
63: *>
64: *> \param[in,out] A
65: *> \verbatim
66: *> A is DOUBLE PRECISION array, dimension (LDA,N)
67: *> On entry, the M-by-N matrix A.
68: *> On exit, the elements on and above the diagonal of the array
69: *> contain the min(M,N)-by-N upper trapezoidal matrix R (R is
70: *> upper triangular if m >= n); the elements below the diagonal,
71: *> with the array TAU, represent the orthogonal matrix Q as a
72: *> product of min(m,n) elementary reflectors (see Further
73: *> Details).
74: *> \endverbatim
75: *>
76: *> \param[in] LDA
77: *> \verbatim
78: *> LDA is INTEGER
79: *> The leading dimension of the array A. LDA >= max(1,M).
80: *> \endverbatim
81: *>
82: *> \param[out] TAU
83: *> \verbatim
84: *> TAU is DOUBLE PRECISION array, dimension (min(M,N))
85: *> The scalar factors of the elementary reflectors (see Further
86: *> Details).
87: *> \endverbatim
88: *>
89: *> \param[out] WORK
90: *> \verbatim
91: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
92: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
93: *> \endverbatim
94: *>
95: *> \param[in] LWORK
96: *> \verbatim
97: *> LWORK is INTEGER
98: *> The dimension of the array WORK. LWORK >= max(1,N).
99: *> For optimum performance LWORK >= N*NB, where NB is
100: *> the optimal blocksize.
101: *>
102: *> If LWORK = -1, then a workspace query is assumed; the routine
103: *> only calculates the optimal size of the WORK array, returns
104: *> this value as the first entry of the WORK array, and no error
105: *> message related to LWORK is issued by XERBLA.
106: *> \endverbatim
107: *>
108: *> \param[out] INFO
109: *> \verbatim
110: *> INFO is INTEGER
111: *> = 0: successful exit
112: *> < 0: if INFO = -i, the i-th argument had an illegal value
113: *> \endverbatim
114: *
115: * Authors:
116: * ========
117: *
1.15 bertrand 118: *> \author Univ. of Tennessee
119: *> \author Univ. of California Berkeley
120: *> \author Univ. of Colorado Denver
121: *> \author NAG Ltd.
1.9 bertrand 122: *
1.18 ! bertrand 123: *> \date November 2019
1.9 bertrand 124: *
125: *> \ingroup doubleGEcomputational
126: *
127: *> \par Further Details:
128: * =====================
129: *>
130: *> \verbatim
131: *>
132: *> The matrix Q is represented as a product of elementary reflectors
133: *>
134: *> Q = H(1) H(2) . . . H(k), where k = min(m,n).
135: *>
136: *> Each H(i) has the form
137: *>
138: *> H(i) = I - tau * v * v**T
139: *>
140: *> where tau is a real scalar, and v is a real vector with
141: *> v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
142: *> and tau in TAU(i).
143: *> \endverbatim
144: *>
145: * =====================================================================
1.1 bertrand 146: SUBROUTINE DGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
147: *
1.18 ! bertrand 148: * -- LAPACK computational routine (version 3.9.0) --
1.1 bertrand 149: * -- LAPACK is a software package provided by Univ. of Tennessee, --
150: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.18 ! bertrand 151: * November 2019
1.1 bertrand 152: *
153: * .. Scalar Arguments ..
154: INTEGER INFO, LDA, LWORK, M, N
155: * ..
156: * .. Array Arguments ..
157: DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
158: * ..
159: *
160: * =====================================================================
161: *
162: * .. Local Scalars ..
163: LOGICAL LQUERY
164: INTEGER I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB,
165: $ NBMIN, NX
166: * ..
167: * .. External Subroutines ..
168: EXTERNAL DGEQR2, DLARFB, DLARFT, XERBLA
169: * ..
170: * .. Intrinsic Functions ..
171: INTRINSIC MAX, MIN
172: * ..
173: * .. External Functions ..
174: INTEGER ILAENV
175: EXTERNAL ILAENV
176: * ..
177: * .. Executable Statements ..
178: *
179: * Test the input arguments
180: *
181: INFO = 0
182: NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
183: LWKOPT = N*NB
184: WORK( 1 ) = LWKOPT
185: LQUERY = ( LWORK.EQ.-1 )
186: IF( M.LT.0 ) THEN
187: INFO = -1
188: ELSE IF( N.LT.0 ) THEN
189: INFO = -2
190: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
191: INFO = -4
192: ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
193: INFO = -7
194: END IF
195: IF( INFO.NE.0 ) THEN
196: CALL XERBLA( 'DGEQRF', -INFO )
197: RETURN
198: ELSE IF( LQUERY ) THEN
199: RETURN
200: END IF
201: *
202: * Quick return if possible
203: *
204: K = MIN( M, N )
205: IF( K.EQ.0 ) THEN
206: WORK( 1 ) = 1
207: RETURN
208: END IF
209: *
210: NBMIN = 2
211: NX = 0
212: IWS = N
213: IF( NB.GT.1 .AND. NB.LT.K ) THEN
214: *
215: * Determine when to cross over from blocked to unblocked code.
216: *
217: NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) )
218: IF( NX.LT.K ) THEN
219: *
220: * Determine if workspace is large enough for blocked code.
221: *
222: LDWORK = N
223: IWS = LDWORK*NB
224: IF( LWORK.LT.IWS ) THEN
225: *
226: * Not enough workspace to use optimal NB: reduce NB and
227: * determine the minimum value of NB.
228: *
229: NB = LWORK / LDWORK
230: NBMIN = MAX( 2, ILAENV( 2, 'DGEQRF', ' ', M, N, -1,
231: $ -1 ) )
232: END IF
233: END IF
234: END IF
235: *
236: IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
237: *
238: * Use blocked code initially
239: *
240: DO 10 I = 1, K - NX, NB
241: IB = MIN( K-I+1, NB )
242: *
243: * Compute the QR factorization of the current block
244: * A(i:m,i:i+ib-1)
245: *
246: CALL DGEQR2( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK,
247: $ IINFO )
248: IF( I+IB.LE.N ) THEN
249: *
250: * Form the triangular factor of the block reflector
251: * H = H(i) H(i+1) . . . H(i+ib-1)
252: *
253: CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB,
254: $ A( I, I ), LDA, TAU( I ), WORK, LDWORK )
255: *
1.8 bertrand 256: * Apply H**T to A(i:m,i+ib:n) from the left
1.1 bertrand 257: *
258: CALL DLARFB( 'Left', 'Transpose', 'Forward',
259: $ 'Columnwise', M-I+1, N-I-IB+1, IB,
260: $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ),
261: $ LDA, WORK( IB+1 ), LDWORK )
262: END IF
263: 10 CONTINUE
264: ELSE
265: I = 1
266: END IF
267: *
268: * Use unblocked code to factor the last or only block.
269: *
270: IF( I.LE.K )
271: $ CALL DGEQR2( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK,
272: $ IINFO )
273: *
274: WORK( 1 ) = IWS
275: RETURN
276: *
277: * End of DGEQRF
278: *
279: END
CVSweb interface <joel.bertrand@systella.fr>