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