LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaorhr_col_getrfnp2.f
Go to the documentation of this file.
1*> \brief \b DLAORHR_COL_GETRFNP2
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLAORHR_GETRF2NP + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* RECURSIVE SUBROUTINE DLAORHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, LDA, M, N
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION A( LDA, * ), D( * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> DLAORHR_COL_GETRFNP2 computes the modified LU factorization without
37*> pivoting of a real 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 DORHR_COL. In DORHR_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*> DLAORHR_COL_GETRFNP2 is called to factorize a block by the blocked
86*> routine DLAORHR_COL_GETRFNP, which uses blocked code calling
87*> Level 3 BLAS to update the submatrix. However, DLAORHR_COL_GETRFNP2
88*> is self-sufficient and can be used without DLAORHR_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 DOUBLE PRECISION 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 DOUBLE PRECISION 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
134*> be only plus or minus one.
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 dlaorhr_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 DOUBLE PRECISION a( lda, * ), d( * )
178* ..
179*
180* =====================================================================
181*
182* .. Parameters ..
183 DOUBLE PRECISION one
184 parameter( one = 1.0d+0 )
185* ..
186* .. Local Scalars ..
187 DOUBLE PRECISION sfmin
188 INTEGER i, iinfo, n1, n2
189* ..
190* .. External Functions ..
191 DOUBLE PRECISION dlamch
192 EXTERNAL dlamch
193* ..
194* .. External Subroutines ..
195 EXTERNAL dgemm, dscal, dtrsm, xerbla
196* ..
197* .. Intrinsic Functions ..
198 INTRINSIC abs, dsign, max, min
199* ..
200* .. Executable Statements ..
201*
202* Test the input parameters
203*
204 info = 0
205 IF( m.LT.0 ) THEN
206 info = -1
207 ELSE IF( n.LT.0 ) THEN
208 info = -2
209 ELSE IF( lda.LT.max( 1, m ) ) THEN
210 info = -4
211 END IF
212 IF( info.NE.0 ) THEN
213 CALL xerbla( 'DLAORHR_COL_GETRFNP2', -info )
214 RETURN
215 END IF
216*
217* Quick return if possible
218*
219 IF( min( m, n ).EQ.0 )
220 $ RETURN
221
222 IF ( m.EQ.1 ) THEN
223*
224* One row case, (also recursion termination case),
225* use unblocked code
226*
227* Transfer the sign
228*
229 d( 1 ) = -dsign( one, a( 1, 1 ) )
230*
231* Construct the row of U
232*
233 a( 1, 1 ) = a( 1, 1 ) - d( 1 )
234*
235 ELSE IF( n.EQ.1 ) THEN
236*
237* One column case, (also recursion termination case),
238* use unblocked code
239*
240* Transfer the sign
241*
242 d( 1 ) = -dsign( one, a( 1, 1 ) )
243*
244* Construct the row of U
245*
246 a( 1, 1 ) = a( 1, 1 ) - d( 1 )
247*
248* Scale the elements 2:M of the column
249*
250* Determine machine safe minimum
251*
252 sfmin = dlamch('S')
253*
254* Construct the subdiagonal elements of L
255*
256 IF( abs( a( 1, 1 ) ) .GE. sfmin ) THEN
257 CALL dscal( m-1, one / a( 1, 1 ), a( 2, 1 ), 1 )
258 ELSE
259 DO i = 2, m
260 a( i, 1 ) = a( i, 1 ) / a( 1, 1 )
261 END DO
262 END IF
263*
264 ELSE
265*
266* Divide the matrix B into four submatrices
267*
268 n1 = min( m, n ) / 2
269 n2 = n-n1
270
271*
272* Factor B11, recursive call
273*
274 CALL dlaorhr_col_getrfnp2( n1, n1, a, lda, d, iinfo )
275*
276* Solve for B21
277*
278 CALL dtrsm( 'R', 'U', 'N', 'N', m-n1, n1, one, a, lda,
279 $ a( n1+1, 1 ), lda )
280*
281* Solve for B12
282*
283 CALL dtrsm( 'L', 'L', 'N', 'U', n1, n2, one, a, lda,
284 $ a( 1, n1+1 ), lda )
285*
286* Update B22, i.e. compute the Schur complement
287* B22 := B22 - B21*B12
288*
289 CALL dgemm( 'N', 'N', m-n1, n2, n1, -one, a( n1+1, 1 ), lda,
290 $ a( 1, n1+1 ), lda, one, a( n1+1, n1+1 ), lda )
291*
292* Factor B22, recursive call
293*
294 CALL dlaorhr_col_getrfnp2( m-n1, n2, a( n1+1, n1+1 ), lda,
295 $ d( n1+1 ), iinfo )
296*
297 END IF
298 RETURN
299*
300* End of DLAORHR_COL_GETRFNP2
301*
302 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
recursive subroutine dlaorhr_col_getrfnp2(m, n, a, lda, d, info)
DLAORHR_COL_GETRFNP2
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181