LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zlaunhr_col_getrfnp2.f
Go to the documentation of this file.
1*> \brief \b ZLAUNHR_COL_GETRFNP2
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* RECURSIVE SUBROUTINE ZLAUNHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, LDA, M, N
25* ..
26* .. Array Arguments ..
27* COMPLEX*16 A( LDA, * ), D( * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> ZLAUNHR_COL_GETRFNP2 computes the modified LU factorization without
37*> pivoting of a complex general M-by-N matrix A. The factorization has
38*> the form:
39*>
40*> A - S = L * U,
41*>
42*> where:
43*> S is a m-by-n diagonal sign matrix with the diagonal D, so that
44*> D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
45*> as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
46*> i-1 steps of Gaussian elimination. This means that the diagonal
47*> element at each step of "modified" Gaussian elimination is at
48*> least one in absolute value (so that division-by-zero not
49*> possible during the division by the diagonal element);
50*>
51*> L is a M-by-N lower triangular matrix with unit diagonal elements
52*> (lower trapezoidal if M > N);
53*>
54*> and U is a M-by-N upper triangular matrix
55*> (upper trapezoidal if M < N).
56*>
57*> This routine is an auxiliary routine used in the Householder
58*> reconstruction routine ZUNHR_COL. In ZUNHR_COL, this routine is
59*> applied to an M-by-N matrix A with orthonormal columns, where each
60*> element is bounded by one in absolute value. With the choice of
61*> the matrix S above, one can show that the diagonal element at each
62*> step of Gaussian elimination is the largest (in absolute value) in
63*> the column on or below the diagonal, so that no pivoting is required
64*> for numerical stability [1].
65*>
66*> For more details on the Householder reconstruction algorithm,
67*> including the modified LU factorization, see [1].
68*>
69*> This is the recursive version of the LU factorization algorithm.
70*> Denote A - S by B. The algorithm divides the matrix B into four
71*> submatrices:
72*>
73*> [ B11 | B12 ] where B11 is n1 by n1,
74*> B = [ -----|----- ] B21 is (m-n1) by n1,
75*> [ B21 | B22 ] B12 is n1 by n2,
76*> B22 is (m-n1) by n2,
77*> with n1 = min(m,n)/2, n2 = n-n1.
78*>
79*>
80*> The subroutine calls itself to factor B11, solves for B21,
81*> solves for B12, updates B22, then calls itself to factor B22.
82*>
83*> For more details on the recursive LU algorithm, see [2].
84*>
85*> ZLAUNHR_COL_GETRFNP2 is called to factorize a block by the blocked
86*> routine ZLAUNHR_COL_GETRFNP, which uses blocked code calling
87*> Level 3 BLAS to update the submatrix. However, ZLAUNHR_COL_GETRFNP2
88*> is self-sufficient and can be used without ZLAUNHR_COL_GETRFNP.
89*>
90*> [1] "Reconstructing Householder vectors from tall-skinny QR",
91*> G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
92*> E. Solomonik, J. Parallel Distrib. Comput.,
93*> vol. 85, pp. 3-31, 2015.
94*>
95*> [2] "Recursion leads to automatic variable blocking for dense linear
96*> algebra algorithms", F. Gustavson, IBM J. of Res. and Dev.,
97*> vol. 41, no. 6, pp. 737-755, 1997.
98*> \endverbatim
99*
100* Arguments:
101* ==========
102*
103*> \param[in] M
104*> \verbatim
105*> M is INTEGER
106*> The number of rows of the matrix A. M >= 0.
107*> \endverbatim
108*>
109*> \param[in] N
110*> \verbatim
111*> N is INTEGER
112*> The number of columns of the matrix A. N >= 0.
113*> \endverbatim
114*>
115*> \param[in,out] A
116*> \verbatim
117*> A is COMPLEX*16 array, dimension (LDA,N)
118*> On entry, the M-by-N matrix to be factored.
119*> On exit, the factors L and U from the factorization
120*> A-S=L*U; the unit diagonal elements of L are not stored.
121*> \endverbatim
122*>
123*> \param[in] LDA
124*> \verbatim
125*> LDA is INTEGER
126*> The leading dimension of the array A. LDA >= max(1,M).
127*> \endverbatim
128*>
129*> \param[out] D
130*> \verbatim
131*> D is COMPLEX*16 array, dimension min(M,N)
132*> The diagonal elements of the diagonal M-by-N sign matrix S,
133*> D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can be
134*> only ( +1.0, 0.0 ) or (-1.0, 0.0 ).
135*> \endverbatim
136*>
137*> \param[out] INFO
138*> \verbatim
139*> INFO is INTEGER
140*> = 0: successful exit
141*> < 0: if INFO = -i, the i-th argument had an illegal value
142*> \endverbatim
143*>
144* Authors:
145* ========
146*
147*> \author Univ. of Tennessee
148*> \author Univ. of California Berkeley
149*> \author Univ. of Colorado Denver
150*> \author NAG Ltd.
151*
152*> \ingroup launhr_col_getrfnp2
153*
154*> \par Contributors:
155* ==================
156*>
157*> \verbatim
158*>
159*> November 2019, Igor Kozachenko,
160*> Computer Science Division,
161*> University of California, Berkeley
162*>
163*> \endverbatim
164*
165* =====================================================================
166 RECURSIVE SUBROUTINE zlaunhr_col_getrfnp2( M, N, A, LDA, D, INFO )
167 IMPLICIT NONE
168*
169* -- LAPACK computational routine --
170* -- LAPACK is a software package provided by Univ. of Tennessee, --
171* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
172*
173* .. Scalar Arguments ..
174 INTEGER info, lda, m, n
175* ..
176* .. Array Arguments ..
177 COMPLEX*16 a( lda, * ), d( * )
178* ..
179*
180* =====================================================================
181*
182* .. Parameters ..
183 DOUBLE PRECISION one
184 parameter( one = 1.0d+0 )
185 COMPLEX*16 cone
186 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
187* ..
188* .. Local Scalars ..
189 DOUBLE PRECISION sfmin
190 INTEGER i, iinfo, n1, n2
191 COMPLEX*16 z
192* ..
193* .. External Functions ..
194 DOUBLE PRECISION dlamch
195 EXTERNAL dlamch
196* ..
197* .. External Subroutines ..
198 EXTERNAL zgemm, zscal, ztrsm, xerbla
199* ..
200* .. Intrinsic Functions ..
201 INTRINSIC abs, dble, dcmplx, dimag, dsign, max, min
202* ..
203* .. Statement Functions ..
204 DOUBLE PRECISION cabs1
205* ..
206* .. Statement Function definitions ..
207 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
208* ..
209* .. Executable Statements ..
210*
211* Test the input parameters
212*
213 info = 0
214 IF( m.LT.0 ) THEN
215 info = -1
216 ELSE IF( n.LT.0 ) THEN
217 info = -2
218 ELSE IF( lda.LT.max( 1, m ) ) THEN
219 info = -4
220 END IF
221 IF( info.NE.0 ) THEN
222 CALL xerbla( 'ZLAUNHR_COL_GETRFNP2', -info )
223 RETURN
224 END IF
225*
226* Quick return if possible
227*
228 IF( min( m, n ).EQ.0 )
229 \$ RETURN
230
231 IF ( m.EQ.1 ) THEN
232*
233* One row case, (also recursion termination case),
234* use unblocked code
235*
236* Transfer the sign
237*
238 d( 1 ) = dcmplx( -dsign( one, dble( a( 1, 1 ) ) ) )
239*
240* Construct the row of U
241*
242 a( 1, 1 ) = a( 1, 1 ) - d( 1 )
243*
244 ELSE IF( n.EQ.1 ) THEN
245*
246* One column case, (also recursion termination case),
247* use unblocked code
248*
249* Transfer the sign
250*
251 d( 1 ) = dcmplx( -dsign( one, dble( a( 1, 1 ) ) ) )
252*
253* Construct the row of U
254*
255 a( 1, 1 ) = a( 1, 1 ) - d( 1 )
256*
257* Scale the elements 2:M of the column
258*
259* Determine machine safe minimum
260*
261 sfmin = dlamch('S')
262*
263* Construct the subdiagonal elements of L
264*
265 IF( cabs1( a( 1, 1 ) ) .GE. sfmin ) THEN
266 CALL zscal( m-1, cone / a( 1, 1 ), a( 2, 1 ), 1 )
267 ELSE
268 DO i = 2, m
269 a( i, 1 ) = a( i, 1 ) / a( 1, 1 )
270 END DO
271 END IF
272*
273 ELSE
274*
275* Divide the matrix B into four submatrices
276*
277 n1 = min( m, n ) / 2
278 n2 = n-n1
279
280*
281* Factor B11, recursive call
282*
283 CALL zlaunhr_col_getrfnp2( n1, n1, a, lda, d, iinfo )
284*
285* Solve for B21
286*
287 CALL ztrsm( 'R', 'U', 'N', 'N', m-n1, n1, cone, a, lda,
288 \$ a( n1+1, 1 ), lda )
289*
290* Solve for B12
291*
292 CALL ztrsm( 'L', 'L', 'N', 'U', n1, n2, cone, a, lda,
293 \$ a( 1, n1+1 ), lda )
294*
295* Update B22, i.e. compute the Schur complement
296* B22 := B22 - B21*B12
297*
298 CALL zgemm( 'N', 'N', m-n1, n2, n1, -cone, a( n1+1, 1 ), lda,
299 \$ a( 1, n1+1 ), lda, cone, a( n1+1, n1+1 ), lda )
300*
301* Factor B22, recursive call
302*
303 CALL zlaunhr_col_getrfnp2( m-n1, n2, a( n1+1, n1+1 ), lda,
304 \$ d( n1+1 ), iinfo )
305*
306 END IF
307 RETURN
308*
309* End of ZLAUNHR_COL_GETRFNP2
310*
311 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
recursive subroutine zlaunhr_col_getrfnp2(m, n, a, lda, d, info)
ZLAUNHR_COL_GETRFNP2
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180