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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtzrzf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack