1*> \brief \b DGEQRT
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGEQRT + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqrt.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqrt.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqrt.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DGEQRT( M, N, NB, A, LDA, T, LDT, WORK, INFO )
22*
23*       .. Scalar Arguments ..
24*       INTEGER INFO, LDA, LDT, M, N, NB
25*       ..
26*       .. Array Arguments ..
27*       DOUBLE PRECISION A( LDA, * ), T( LDT, * ), WORK( * )
28*       ..
29*
30*
31*> \par Purpose:
32*  =============
33*>
34*> \verbatim
35*>
36*> DGEQRT computes a blocked QR factorization of a real M-by-N matrix A
37*> using the compact WY representation of Q.
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] NB
56*> \verbatim
57*>          NB is INTEGER
58*>          The block size to be used in the blocked QR.  MIN(M,N) >= NB >= 1.
59*> \endverbatim
60*>
61*> \param[in,out] A
62*> \verbatim
63*>          A is DOUBLE PRECISION array, dimension (LDA,N)
64*>          On entry, the M-by-N matrix A.
65*>          On exit, the elements on and above the diagonal of the array
66*>          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
67*>          upper triangular if M >= N); the elements below the diagonal
68*>          are the columns of V.
69*> \endverbatim
70*>
71*> \param[in] LDA
72*> \verbatim
73*>          LDA is INTEGER
74*>          The leading dimension of the array A.  LDA >= max(1,M).
75*> \endverbatim
76*>
77*> \param[out] T
78*> \verbatim
79*>          T is DOUBLE PRECISION array, dimension (LDT,MIN(M,N))
80*>          The upper triangular block reflectors stored in compact form
81*>          as a sequence of upper triangular blocks.  See below
82*>          for further details.
83*> \endverbatim
84*>
85*> \param[in] LDT
86*> \verbatim
87*>          LDT is INTEGER
88*>          The leading dimension of the array T.  LDT >= NB.
89*> \endverbatim
90*>
91*> \param[out] WORK
92*> \verbatim
93*>          WORK is DOUBLE PRECISION array, dimension (NB*N)
94*> \endverbatim
95*>
96*> \param[out] INFO
97*> \verbatim
98*>          INFO is INTEGER
99*>          = 0:  successful exit
100*>          < 0:  if INFO = -i, the i-th argument had an illegal value
101*> \endverbatim
102*
103*  Authors:
104*  ========
105*
106*> \author Univ. of Tennessee
107*> \author Univ. of California Berkeley
108*> \author Univ. of Colorado Denver
109*> \author NAG Ltd.
110*
111*> \ingroup doubleGEcomputational
112*
113*> \par Further Details:
114*  =====================
115*>
116*> \verbatim
117*>
118*>  The matrix V stores the elementary reflectors H(i) in the i-th column
119*>  below the diagonal. For example, if M=5 and N=3, the matrix V is
120*>
121*>               V = (  1       )
122*>                   ( v1  1    )
123*>                   ( v1 v2  1 )
124*>                   ( v1 v2 v3 )
125*>                   ( v1 v2 v3 )
126*>
127*>  where the vi's represent the vectors which define H(i), which are returned
128*>  in the matrix A.  The 1's along the diagonal of V are not stored in A.
129*>
130*>  Let K=MIN(M,N).  The number of blocks is B = ceiling(K/NB), where each
131*>  block is of order NB except for the last block, which is of order
132*>  IB = K - (B-1)*NB.  For each of the B blocks, a upper triangular block
133*>  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB
134*>  for the last block) T's are stored in the NB-by-K matrix T as
135*>
136*>               T = (T1 T2 ... TB).
137*> \endverbatim
138*>
139*  =====================================================================
140      SUBROUTINE DGEQRT( M, N, NB, A, LDA, T, LDT, WORK, INFO )
141*
142*  -- LAPACK computational routine --
143*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
144*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145*
146*     .. Scalar Arguments ..
147      INTEGER INFO, LDA, LDT, M, N, NB
148*     ..
149*     .. Array Arguments ..
150      DOUBLE PRECISION A( LDA, * ), T( LDT, * ), WORK( * )
151*     ..
152*
153* =====================================================================
154*
155*     ..
156*     .. Local Scalars ..
157      INTEGER    I, IB, IINFO, K
158      LOGICAL    USE_RECURSIVE_QR
159      PARAMETER( USE_RECURSIVE_QR=.TRUE. )
160*     ..
161*     .. External Subroutines ..
162      EXTERNAL   DGEQRT2, DGEQRT3, DLARFB, XERBLA
163*     ..
164*     .. Executable Statements ..
165*
166*     Test the input arguments
167*
168      INFO = 0
169      IF( M.LT.0 ) THEN
170         INFO = -1
171      ELSE IF( N.LT.0 ) THEN
172         INFO = -2
173      ELSE IF( NB.LT.1 .OR. ( NB.GT.MIN(M,N) .AND. MIN(M,N).GT.0 ) )THEN
174         INFO = -3
175      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
176         INFO = -5
177      ELSE IF( LDT.LT.NB ) THEN
178         INFO = -7
179      END IF
180      IF( INFO.NE.0 ) THEN
181         CALL XERBLA( 'DGEQRT', -INFO )
182         RETURN
183      END IF
184*
185*     Quick return if possible
186*
187      K = MIN( M, N )
188      IF( K.EQ.0 ) RETURN
189*
190*     Blocked loop of length K
191*
192      DO I = 1, K,  NB
193         IB = MIN( K-I+1, NB )
194*
195*     Compute the QR factorization of the current block A(I:M,I:I+IB-1)
196*
197         IF( USE_RECURSIVE_QR ) THEN
198            CALL DGEQRT3( M-I+1, IB, A(I,I), LDA, T(1,I), LDT, IINFO )
199         ELSE
200            CALL DGEQRT2( M-I+1, IB, A(I,I), LDA, T(1,I), LDT, IINFO )
201         END IF
202         IF( I+IB.LE.N ) THEN
203*
204*     Update by applying H**T to A(I:M,I+IB:N) from the left
205*
206            CALL DLARFB( 'L', 'T', 'F', 'C', M-I+1, N-I-IB+1, IB,
207     $                   A( I, I ), LDA, T( 1, I ), LDT,
208     $                   A( I, I+IB ), LDA, WORK , N-I-IB+1 )
209         END IF
210      END DO
211      RETURN
212*
213*     End of DGEQRT
214*
215      END
216