LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date April 2012
120 *
121 *> \ingroup doubleOTHERcomputational
122 *
123 *> \par Contributors:
124 * ==================
125 *>
126 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
127 *
128 *> \par Further Details:
129 * =====================
130 *>
131 *> \verbatim
132 *>
133 *> The N-by-N matrix Z can be computed by
134 *>
135 *> Z = Z(1)*Z(2)* ... *Z(M)
136 *>
137 *> where each N-by-N Z(k) is given by
138 *>
139 *> Z(k) = I - tau(k)*v(k)*v(k)**T
140 *>
141 *> with v(k) is the kth row vector of the M-by-N matrix
142 *>
143 *> V = ( I A(:,M+1:N) )
144 *>
145 *> I is the M-by-M identity matrix, A(:,M+1:N)
146 *> is the output stored in A on exit from DTZRZF,
147 *> and tau(k) is the kth element of the array TAU.
148 *>
149 *> \endverbatim
150 *>
151 * =====================================================================
152  SUBROUTINE dtzrzf( M, N, A, LDA, TAU, WORK, LWORK, INFO )
153 *
154 * -- LAPACK computational routine (version 3.4.1) --
155 * -- LAPACK is a software package provided by Univ. of Tennessee, --
156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157 * April 2012
158 *
159 * .. Scalar Arguments ..
160  INTEGER info, lda, lwork, m, n
161 * ..
162 * .. Array Arguments ..
163  DOUBLE PRECISION a( lda, * ), tau( * ), work( * )
164 * ..
165 *
166 * =====================================================================
167 *
168 * .. Parameters ..
169  DOUBLE PRECISION zero
170  parameter( zero = 0.0d+0 )
171 * ..
172 * .. Local Scalars ..
173  LOGICAL lquery
174  INTEGER i, ib, iws, ki, kk, ldwork, lwkmin, lwkopt,
175  $ m1, mu, nb, nbmin, nx
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL xerbla, dlarzb, dlarzt, dlatrz
179 * ..
180 * .. Intrinsic Functions ..
181  INTRINSIC max, min
182 * ..
183 * .. External Functions ..
184  INTEGER ilaenv
185  EXTERNAL ilaenv
186 * ..
187 * .. Executable Statements ..
188 *
189 * Test the input arguments
190 *
191  info = 0
192  lquery = ( lwork.EQ.-1 )
193  IF( m.LT.0 ) THEN
194  info = -1
195  ELSE IF( n.LT.m ) THEN
196  info = -2
197  ELSE IF( lda.LT.max( 1, m ) ) THEN
198  info = -4
199  END IF
200 *
201  IF( info.EQ.0 ) THEN
202  IF( m.EQ.0 .OR. m.EQ.n ) THEN
203  lwkopt = 1
204  lwkmin = 1
205  ELSE
206 *
207 * Determine the block size.
208 *
209  nb = ilaenv( 1, 'DGERQF', ' ', m, n, -1, -1 )
210  lwkopt = m*nb
211  lwkmin = max( 1, m )
212  END IF
213  work( 1 ) = lwkopt
214 *
215  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
216  info = -7
217  END IF
218  END IF
219 *
220  IF( info.NE.0 ) THEN
221  CALL xerbla( 'DTZRZF', -info )
222  return
223  ELSE IF( lquery ) THEN
224  return
225  END IF
226 *
227 * Quick return if possible
228 *
229  IF( m.EQ.0 ) THEN
230  return
231  ELSE IF( m.EQ.n ) THEN
232  DO 10 i = 1, n
233  tau( i ) = zero
234  10 continue
235  return
236  END IF
237 *
238  nbmin = 2
239  nx = 1
240  iws = m
241  IF( nb.GT.1 .AND. nb.LT.m ) THEN
242 *
243 * Determine when to cross over from blocked to unblocked code.
244 *
245  nx = max( 0, ilaenv( 3, 'DGERQF', ' ', m, n, -1, -1 ) )
246  IF( nx.LT.m ) THEN
247 *
248 * Determine if workspace is large enough for blocked code.
249 *
250  ldwork = m
251  iws = ldwork*nb
252  IF( lwork.LT.iws ) THEN
253 *
254 * Not enough workspace to use optimal NB: reduce NB and
255 * determine the minimum value of NB.
256 *
257  nb = lwork / ldwork
258  nbmin = max( 2, ilaenv( 2, 'DGERQF', ' ', m, n, -1,
259  $ -1 ) )
260  END IF
261  END IF
262  END IF
263 *
264  IF( nb.GE.nbmin .AND. nb.LT.m .AND. nx.LT.m ) THEN
265 *
266 * Use blocked code initially.
267 * The last kk rows are handled by the block method.
268 *
269  m1 = min( m+1, n )
270  ki = ( ( m-nx-1 ) / nb )*nb
271  kk = min( m, ki+nb )
272 *
273  DO 20 i = m - kk + ki + 1, m - kk + 1, -nb
274  ib = min( m-i+1, nb )
275 *
276 * Compute the TZ factorization of the current block
277 * A(i:i+ib-1,i:n)
278 *
279  CALL dlatrz( ib, n-i+1, n-m, a( i, i ), lda, tau( i ),
280  $ work )
281  IF( i.GT.1 ) THEN
282 *
283 * Form the triangular factor of the block reflector
284 * H = H(i+ib-1) . . . H(i+1) H(i)
285 *
286  CALL dlarzt( 'Backward', 'Rowwise', n-m, ib, a( i, m1 ),
287  $ lda, tau( i ), work, ldwork )
288 *
289 * Apply H to A(1:i-1,i:n) from the right
290 *
291  CALL dlarzb( 'Right', 'No transpose', 'Backward',
292  $ 'Rowwise', i-1, n-i+1, ib, n-m, a( i, m1 ),
293  $ lda, work, ldwork, a( 1, i ), lda,
294  $ work( ib+1 ), ldwork )
295  END IF
296  20 continue
297  mu = i + nb - 1
298  ELSE
299  mu = m
300  END IF
301 *
302 * Use unblocked code to factor the last or only block
303 *
304  IF( mu.GT.0 )
305  $ CALL dlatrz( mu, n, n-m, a, lda, tau, work )
306 *
307  work( 1 ) = lwkopt
308 *
309  return
310 *
311 * End of DTZRZF
312 *
313  END