LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dtzrzf.f
Go to the documentation of this file.
1 *> \brief \b DTZRZF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DTZRZF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtzrzf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtzrzf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtzrzf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DTZRZF( 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 *> DTZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A
37 *> to upper triangular form by means of orthogonal transformations.
38 *>
39 *> The upper trapezoidal matrix A is factored as
40 *>
41 *> A = ( R 0 ) * Z,
42 *>
43 *> where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
44 *> triangular matrix.
45 *> \endverbatim
46 *
47 * Arguments:
48 * ==========
49 *
50 *> \param[in] M
51 *> \verbatim
52 *> M is INTEGER
53 *> The number of rows of the matrix A. M >= 0.
54 *> \endverbatim
55 *>
56 *> \param[in] N
57 *> \verbatim
58 *> N is INTEGER
59 *> The number of columns of the matrix A. N >= M.
60 *> \endverbatim
61 *>
62 *> \param[in,out] A
63 *> \verbatim
64 *> A is DOUBLE PRECISION array, dimension (LDA,N)
65 *> On entry, the leading M-by-N upper trapezoidal part of the
66 *> array A must contain the matrix to be factorized.
67 *> On exit, the leading M-by-M upper triangular part of A
68 *> contains the upper triangular matrix R, and elements M+1 to
69 *> N of the first M rows of A, with the array TAU, represent the
70 *> orthogonal matrix Z as a product of M elementary reflectors.
71 *> \endverbatim
72 *>
73 *> \param[in] LDA
74 *> \verbatim
75 *> LDA is INTEGER
76 *> The leading dimension of the array A. LDA >= max(1,M).
77 *> \endverbatim
78 *>
79 *> \param[out] TAU
80 *> \verbatim
81 *> TAU is DOUBLE PRECISION array, dimension (M)
82 *> The scalar factors of the elementary reflectors.
83 *> \endverbatim
84 *>
85 *> \param[out] WORK
86 *> \verbatim
87 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
88 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
89 *> \endverbatim
90 *>
91 *> \param[in] LWORK
92 *> \verbatim
93 *> LWORK is INTEGER
94 *> The dimension of the array WORK. LWORK >= max(1,M).
95 *> For optimum performance LWORK >= M*NB, where NB is
96 *> the optimal blocksize.
97 *>
98 *> If LWORK = -1, then a workspace query is assumed; the routine
99 *> only calculates the optimal size of the WORK array, returns
100 *> this value as the first entry of the WORK array, and no error
101 *> message related to LWORK is issued by XERBLA.
102 *> \endverbatim
103 *>
104 *> \param[out] INFO
105 *> \verbatim
106 *> INFO is INTEGER
107 *> = 0: successful exit
108 *> < 0: if INFO = -i, the i-th argument had an illegal value
109 *> \endverbatim
110 *
111 * Authors:
112 * ========
113 *
114 *> \author Univ. of Tennessee
115 *> \author Univ. of California Berkeley
116 *> \author Univ. of Colorado Denver
117 *> \author NAG Ltd.
118 *
119 *> \ingroup doubleOTHERcomputational
120 *
121 *> \par Contributors:
122 * ==================
123 *>
124 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
125 *
126 *> \par Further Details:
127 * =====================
128 *>
129 *> \verbatim
130 *>
131 *> The N-by-N matrix Z can be computed by
132 *>
133 *> Z = Z(1)*Z(2)* ... *Z(M)
134 *>
135 *> where each N-by-N Z(k) is given by
136 *>
137 *> Z(k) = I - tau(k)*v(k)*v(k)**T
138 *>
139 *> with v(k) is the kth row vector of the M-by-N matrix
140 *>
141 *> V = ( I A(:,M+1:N) )
142 *>
143 *> I is the M-by-M identity matrix, A(:,M+1:N)
144 *> is the output stored in A on exit from DTZRZF,
145 *> and tau(k) is the kth element of the array TAU.
146 *>
147 *> \endverbatim
148 *>
149 * =====================================================================
150  SUBROUTINE dtzrzf( M, N, A, LDA, TAU, WORK, LWORK, INFO )
151 *
152 * -- LAPACK computational routine --
153 * -- LAPACK is a software package provided by Univ. of Tennessee, --
154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155 *
156 * .. Scalar Arguments ..
157  INTEGER INFO, LDA, LWORK, M, N
158 * ..
159 * .. Array Arguments ..
160  DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
161 * ..
162 *
163 * =====================================================================
164 *
165 * .. Parameters ..
166  DOUBLE PRECISION ZERO
167  parameter( zero = 0.0d+0 )
168 * ..
169 * .. Local Scalars ..
170  LOGICAL LQUERY
171  INTEGER I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT,
172  $ M1, MU, NB, NBMIN, NX
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL xerbla, dlarzb, dlarzt, dlatrz
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC max, min
179 * ..
180 * .. External Functions ..
181  INTEGER ILAENV
182  EXTERNAL ilaenv
183 * ..
184 * .. Executable Statements ..
185 *
186 * Test the input arguments
187 *
188  info = 0
189  lquery = ( lwork.EQ.-1 )
190  IF( m.LT.0 ) THEN
191  info = -1
192  ELSE IF( n.LT.m ) THEN
193  info = -2
194  ELSE IF( lda.LT.max( 1, m ) ) THEN
195  info = -4
196  END IF
197 *
198  IF( info.EQ.0 ) THEN
199  IF( m.EQ.0 .OR. m.EQ.n ) THEN
200  lwkopt = 1
201  lwkmin = 1
202  ELSE
203 *
204 * Determine the block size.
205 *
206  nb = ilaenv( 1, 'DGERQF', ' ', m, n, -1, -1 )
207  lwkopt = m*nb
208  lwkmin = max( 1, m )
209  END IF
210  work( 1 ) = lwkopt
211 *
212  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
213  info = -7
214  END IF
215  END IF
216 *
217  IF( info.NE.0 ) THEN
218  CALL xerbla( 'DTZRZF', -info )
219  RETURN
220  ELSE IF( lquery ) THEN
221  RETURN
222  END IF
223 *
224 * Quick return if possible
225 *
226  IF( m.EQ.0 ) THEN
227  RETURN
228  ELSE IF( m.EQ.n ) THEN
229  DO 10 i = 1, n
230  tau( i ) = zero
231  10 CONTINUE
232  RETURN
233  END IF
234 *
235  nbmin = 2
236  nx = 1
237  iws = m
238  IF( nb.GT.1 .AND. nb.LT.m ) THEN
239 *
240 * Determine when to cross over from blocked to unblocked code.
241 *
242  nx = max( 0, ilaenv( 3, 'DGERQF', ' ', m, n, -1, -1 ) )
243  IF( nx.LT.m ) THEN
244 *
245 * Determine if workspace is large enough for blocked code.
246 *
247  ldwork = m
248  iws = ldwork*nb
249  IF( lwork.LT.iws ) THEN
250 *
251 * Not enough workspace to use optimal NB: reduce NB and
252 * determine the minimum value of NB.
253 *
254  nb = lwork / ldwork
255  nbmin = max( 2, ilaenv( 2, 'DGERQF', ' ', m, n, -1,
256  $ -1 ) )
257  END IF
258  END IF
259  END IF
260 *
261  IF( nb.GE.nbmin .AND. nb.LT.m .AND. nx.LT.m ) THEN
262 *
263 * Use blocked code initially.
264 * The last kk rows are handled by the block method.
265 *
266  m1 = min( m+1, n )
267  ki = ( ( m-nx-1 ) / nb )*nb
268  kk = min( m, ki+nb )
269 *
270  DO 20 i = m - kk + ki + 1, m - kk + 1, -nb
271  ib = min( m-i+1, nb )
272 *
273 * Compute the TZ factorization of the current block
274 * A(i:i+ib-1,i:n)
275 *
276  CALL dlatrz( ib, n-i+1, n-m, a( i, i ), lda, tau( i ),
277  $ work )
278  IF( i.GT.1 ) THEN
279 *
280 * Form the triangular factor of the block reflector
281 * H = H(i+ib-1) . . . H(i+1) H(i)
282 *
283  CALL dlarzt( 'Backward', 'Rowwise', n-m, ib, a( i, m1 ),
284  $ lda, tau( i ), work, ldwork )
285 *
286 * Apply H to A(1:i-1,i:n) from the right
287 *
288  CALL dlarzb( 'Right', 'No transpose', 'Backward',
289  $ 'Rowwise', i-1, n-i+1, ib, n-m, a( i, m1 ),
290  $ lda, work, ldwork, a( 1, i ), lda,
291  $ work( ib+1 ), ldwork )
292  END IF
293  20 CONTINUE
294  mu = i + nb - 1
295  ELSE
296  mu = m
297  END IF
298 *
299 * Use unblocked code to factor the last or only block
300 *
301  IF( mu.GT.0 )
302  $ CALL dlatrz( mu, n, n-m, a, lda, tau, work )
303 *
304  work( 1 ) = lwkopt
305 *
306  RETURN
307 *
308 * End of DTZRZF
309 *
310  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlarzt(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
DLARZT forms the triangular factor T of a block reflector H = I - vtvH.
Definition: dlarzt.f:185
subroutine dtzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DTZRZF
Definition: dtzrzf.f:151
subroutine dlatrz(M, N, L, A, LDA, TAU, WORK)
DLATRZ factors an upper trapezoidal matrix by means of orthogonal transformations.
Definition: dlatrz.f:140
subroutine dlarzb(SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, C, LDC, WORK, LDWORK)
DLARZB applies a block reflector or its transpose to a general matrix.
Definition: dlarzb.f:183