LAPACK  3.10.1 LAPACK: Linear Algebra PACKage
dgeqr.f
Go to the documentation of this file.
1 *> \brief \b DGEQR
2 *
3 * Definition:
4 * ===========
5 *
6 * SUBROUTINE DGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK,
7 * INFO )
8 *
9 * .. Scalar Arguments ..
10 * INTEGER INFO, LDA, M, N, TSIZE, LWORK
11 * ..
12 * .. Array Arguments ..
13 * DOUBLE PRECISION A( LDA, * ), T( * ), WORK( * )
14 * ..
15 *
16 *
17 *> \par Purpose:
18 * =============
19 *>
20 *> \verbatim
21 *>
22 *> DGEQR computes a QR factorization of a real M-by-N matrix A:
23 *>
24 *> A = Q * ( R ),
25 *> ( 0 )
26 *>
27 *> where:
28 *>
29 *> Q is a M-by-M orthogonal matrix;
30 *> R is an upper-triangular N-by-N matrix;
31 *> 0 is a (M-N)-by-N zero matrix, if M > N.
32 *>
33 *> \endverbatim
34 *
35 * Arguments:
36 * ==========
37 *
38 *> \param[in] M
39 *> \verbatim
40 *> M is INTEGER
41 *> The number of rows of the matrix A. M >= 0.
42 *> \endverbatim
43 *>
44 *> \param[in] N
45 *> \verbatim
46 *> N is INTEGER
47 *> The number of columns of the matrix A. N >= 0.
48 *> \endverbatim
49 *>
50 *> \param[in,out] A
51 *> \verbatim
52 *> A is DOUBLE PRECISION array, dimension (LDA,N)
53 *> On entry, the M-by-N matrix A.
54 *> On exit, the elements on and above the diagonal of the array
55 *> contain the min(M,N)-by-N upper trapezoidal matrix R
56 *> (R is upper triangular if M >= N);
57 *> the elements below the diagonal are used to store part of the
58 *> data structure to represent Q.
59 *> \endverbatim
60 *>
61 *> \param[in] LDA
62 *> \verbatim
63 *> LDA is INTEGER
64 *> The leading dimension of the array A. LDA >= max(1,M).
65 *> \endverbatim
66 *>
67 *> \param[out] T
68 *> \verbatim
69 *> T is DOUBLE PRECISION array, dimension (MAX(5,TSIZE))
70 *> On exit, if INFO = 0, T(1) returns optimal (or either minimal
71 *> or optimal, if query is assumed) TSIZE. See TSIZE for details.
72 *> Remaining T contains part of the data structure used to represent Q.
73 *> If one wants to apply or construct Q, then one needs to keep T
74 *> (in addition to A) and pass it to further subroutines.
75 *> \endverbatim
76 *>
77 *> \param[in] TSIZE
78 *> \verbatim
79 *> TSIZE is INTEGER
80 *> If TSIZE >= 5, the dimension of the array T.
81 *> If TSIZE = -1 or -2, then a workspace query is assumed. The routine
82 *> only calculates the sizes of the T and WORK arrays, returns these
83 *> values as the first entries of the T and WORK arrays, and no error
84 *> message related to T or WORK is issued by XERBLA.
85 *> If TSIZE = -1, the routine calculates optimal size of T for the
86 *> optimum performance and returns this value in T(1).
87 *> If TSIZE = -2, the routine calculates minimal size of T and
88 *> returns this value in T(1).
89 *> \endverbatim
90 *>
91 *> \param[out] WORK
92 *> \verbatim
93 *> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
94 *> On exit, if INFO = 0, WORK(1) contains optimal (or either minimal
95 *> or optimal, if query was assumed) LWORK.
96 *> See LWORK for details.
97 *> \endverbatim
98 *>
99 *> \param[in] LWORK
100 *> \verbatim
101 *> LWORK is INTEGER
102 *> The dimension of the array WORK.
103 *> If LWORK = -1 or -2, then a workspace query is assumed. The routine
104 *> only calculates the sizes of the T and WORK arrays, returns these
105 *> values as the first entries of the T and WORK arrays, and no error
106 *> message related to T or WORK is issued by XERBLA.
107 *> If LWORK = -1, the routine calculates optimal size of WORK for the
108 *> optimal performance and returns this value in WORK(1).
109 *> If LWORK = -2, the routine calculates minimal size of WORK and
110 *> returns this value in WORK(1).
111 *> \endverbatim
112 *>
113 *> \param[out] INFO
114 *> \verbatim
115 *> INFO is INTEGER
116 *> = 0: successful exit
117 *> < 0: if INFO = -i, the i-th argument had an illegal value
118 *> \endverbatim
119 *
120 * Authors:
121 * ========
122 *
123 *> \author Univ. of Tennessee
124 *> \author Univ. of California Berkeley
125 *> \author Univ. of Colorado Denver
126 *> \author NAG Ltd.
127 *
128 *> \par Further Details
129 * ====================
130 *>
131 *> \verbatim
132 *>
133 *> The goal of the interface is to give maximum freedom to the developers for
134 *> creating any QR factorization algorithm they wish. The triangular
135 *> (trapezoidal) R has to be stored in the upper part of A. The lower part of A
136 *> and the array T can be used to store any relevant information for applying or
137 *> constructing the Q factor. The WORK array can safely be discarded after exit.
138 *>
139 *> Caution: One should not expect the sizes of T and WORK to be the same from one
140 *> LAPACK implementation to the other, or even from one execution to the other.
141 *> A workspace query (for T and WORK) is needed at each execution. However,
142 *> for a given execution, the size of T and WORK are fixed and will not change
143 *> from one query to the next.
144 *>
145 *> \endverbatim
146 *>
147 *> \par Further Details particular to this LAPACK implementation:
148 * ==============================================================
149 *>
150 *> \verbatim
151 *>
152 *> These details are particular for this LAPACK implementation. Users should not
153 *> take them for granted. These details may change in the future, and are not likely
154 *> true for another LAPACK implementation. These details are relevant if one wants
155 *> to try to understand the code. They are not part of the interface.
156 *>
157 *> In this version,
158 *>
159 *> T(2): row block size (MB)
160 *> T(3): column block size (NB)
161 *> T(6:TSIZE): data structure needed for Q, computed by
162 *> DLATSQR or DGEQRT
163 *>
164 *> Depending on the matrix dimensions M and N, and row and column
165 *> block sizes MB and NB returned by ILAENV, DGEQR will use either
166 *> DLATSQR (if the matrix is tall-and-skinny) or DGEQRT to compute
167 *> the QR factorization.
168 *>
169 *> \endverbatim
170 *>
171 * =====================================================================
172  SUBROUTINE dgeqr( M, N, A, LDA, T, TSIZE, WORK, LWORK,
173  \$ INFO )
174 *
175 * -- LAPACK computational routine --
176 * -- LAPACK is a software package provided by Univ. of Tennessee, --
177 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
178 *
179 * .. Scalar Arguments ..
180  INTEGER INFO, LDA, M, N, TSIZE, LWORK
181 * ..
182 * .. Array Arguments ..
183  DOUBLE PRECISION A( LDA, * ), T( * ), WORK( * )
184 * ..
185 *
186 * =====================================================================
187 *
188 * ..
189 * .. Local Scalars ..
190  LOGICAL LQUERY, LMINWS, MINT, MINW
191  INTEGER MB, NB, MINTSZ, NBLCKS
192 * ..
193 * .. External Functions ..
194  LOGICAL LSAME
195  EXTERNAL lsame
196 * ..
197 * .. External Subroutines ..
198  EXTERNAL dlatsqr, dgeqrt, xerbla
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC max, min, mod
202 * ..
203 * .. External Functions ..
204  INTEGER ILAENV
205  EXTERNAL ilaenv
206 * ..
207 * .. Executable Statements ..
208 *
209 * Test the input arguments
210 *
211  info = 0
212 *
213  lquery = ( tsize.EQ.-1 .OR. tsize.EQ.-2 .OR.
214  \$ lwork.EQ.-1 .OR. lwork.EQ.-2 )
215 *
216  mint = .false.
217  minw = .false.
218  IF( tsize.EQ.-2 .OR. lwork.EQ.-2 ) THEN
219  IF( tsize.NE.-1 ) mint = .true.
220  IF( lwork.NE.-1 ) minw = .true.
221  END IF
222 *
223 * Determine the block size
224 *
225  IF( min( m, n ).GT.0 ) THEN
226  mb = ilaenv( 1, 'DGEQR ', ' ', m, n, 1, -1 )
227  nb = ilaenv( 1, 'DGEQR ', ' ', m, n, 2, -1 )
228  ELSE
229  mb = m
230  nb = 1
231  END IF
232  IF( mb.GT.m .OR. mb.LE.n ) mb = m
233  IF( nb.GT.min( m, n ) .OR. nb.LT.1 ) nb = 1
234  mintsz = n + 5
235  IF( mb.GT.n .AND. m.GT.n ) THEN
236  IF( mod( m - n, mb - n ).EQ.0 ) THEN
237  nblcks = ( m - n ) / ( mb - n )
238  ELSE
239  nblcks = ( m - n ) / ( mb - n ) + 1
240  END IF
241  ELSE
242  nblcks = 1
243  END IF
244 *
245 * Determine if the workspace size satisfies minimal size
246 *
247  lminws = .false.
248  IF( ( tsize.LT.max( 1, nb*n*nblcks + 5 ) .OR. lwork.LT.nb*n )
249  \$ .AND. ( lwork.GE.n ) .AND. ( tsize.GE.mintsz )
250  \$ .AND. ( .NOT.lquery ) ) THEN
251  IF( tsize.LT.max( 1, nb*n*nblcks + 5 ) ) THEN
252  lminws = .true.
253  nb = 1
254  mb = m
255  END IF
256  IF( lwork.LT.nb*n ) THEN
257  lminws = .true.
258  nb = 1
259  END IF
260  END IF
261 *
262  IF( m.LT.0 ) THEN
263  info = -1
264  ELSE IF( n.LT.0 ) THEN
265  info = -2
266  ELSE IF( lda.LT.max( 1, m ) ) THEN
267  info = -4
268  ELSE IF( tsize.LT.max( 1, nb*n*nblcks + 5 )
269  \$ .AND. ( .NOT.lquery ) .AND. ( .NOT.lminws ) ) THEN
270  info = -6
271  ELSE IF( ( lwork.LT.max( 1, n*nb ) ) .AND. ( .NOT.lquery )
272  \$ .AND. ( .NOT.lminws ) ) THEN
273  info = -8
274  END IF
275 *
276  IF( info.EQ.0 ) THEN
277  IF( mint ) THEN
278  t( 1 ) = mintsz
279  ELSE
280  t( 1 ) = nb*n*nblcks + 5
281  END IF
282  t( 2 ) = mb
283  t( 3 ) = nb
284  IF( minw ) THEN
285  work( 1 ) = max( 1, n )
286  ELSE
287  work( 1 ) = max( 1, nb*n )
288  END IF
289  END IF
290  IF( info.NE.0 ) THEN
291  CALL xerbla( 'DGEQR', -info )
292  RETURN
293  ELSE IF( lquery ) THEN
294  RETURN
295  END IF
296 *
297 * Quick return if possible
298 *
299  IF( min( m, n ).EQ.0 ) THEN
300  RETURN
301  END IF
302 *
303 * The QR Decomposition
304 *
305  IF( ( m.LE.n ) .OR. ( mb.LE.n ) .OR. ( mb.GE.m ) ) THEN
306  CALL dgeqrt( m, n, nb, a, lda, t( 6 ), nb, work, info )
307  ELSE
308  CALL dlatsqr( m, n, mb, nb, a, lda, t( 6 ), nb, work,
309  \$ lwork, info )
310  END IF
311 *
312  work( 1 ) = max( 1, nb*n )
313 *
314  RETURN
315 *
316 * End of DGEQR
317 *
318  END
subroutine dgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
DGEQR
Definition: dgeqr.f:174
subroutine dlatsqr(M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK, INFO)
DLATSQR
Definition: dlatsqr.f:166
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dgeqrt(M, N, NB, A, LDA, T, LDT, WORK, INFO)
DGEQRT
Definition: dgeqrt.f:141