LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlaorhr_col_getrfnp2()

recursive subroutine dlaorhr_col_getrfnp2 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  D,
integer  INFO 
)

DLAORHR_COL_GETRFNP2

Download DLAORHR_GETRF2NP + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLAORHR_COL_GETRFNP2 computes the modified LU factorization without
 pivoting of a real general M-by-N matrix A. The factorization has
 the form:

     A - S = L * U,

 where:
    S is a m-by-n diagonal sign matrix with the diagonal D, so that
    D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
    as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
    i-1 steps of Gaussian elimination. This means that the diagonal
    element at each step of "modified" Gaussian elimination is at
    least one in absolute value (so that division-by-zero not
    possible during the division by the diagonal element);

    L is a M-by-N lower triangular matrix with unit diagonal elements
    (lower trapezoidal if M > N);

    and U is a M-by-N upper triangular matrix
    (upper trapezoidal if M < N).

 This routine is an auxiliary routine used in the Householder
 reconstruction routine DORHR_COL. In DORHR_COL, this routine is
 applied to an M-by-N matrix A with orthonormal columns, where each
 element is bounded by one in absolute value. With the choice of
 the matrix S above, one can show that the diagonal element at each
 step of Gaussian elimination is the largest (in absolute value) in
 the column on or below the diagonal, so that no pivoting is required
 for numerical stability [1].

 For more details on the Householder reconstruction algorithm,
 including the modified LU factorization, see [1].

 This is the recursive version of the LU factorization algorithm.
 Denote A - S by B. The algorithm divides the matrix B into four
 submatrices:

        [  B11 | B12  ]  where B11 is n1 by n1,
    B = [ -----|----- ]        B21 is (m-n1) by n1,
        [  B21 | B22  ]        B12 is n1 by n2,
                               B22 is (m-n1) by n2,
                               with n1 = min(m,n)/2, n2 = n-n1.


 The subroutine calls itself to factor B11, solves for B21,
 solves for B12, updates B22, then calls itself to factor B22.

 For more details on the recursive LU algorithm, see [2].

 DLAORHR_COL_GETRFNP2 is called to factorize a block by the blocked
 routine DLAORHR_COL_GETRFNP, which uses blocked code calling
 Level 3 BLAS to update the submatrix. However, DLAORHR_COL_GETRFNP2
 is self-sufficient and can be used without DLAORHR_COL_GETRFNP.

 [1] "Reconstructing Householder vectors from tall-skinny QR",
     G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
     E. Solomonik, J. Parallel Distrib. Comput.,
     vol. 85, pp. 3-31, 2015.

 [2] "Recursion leads to automatic variable blocking for dense linear
     algebra algorithms", F. Gustavson, IBM J. of Res. and Dev.,
     vol. 41, no. 6, pp. 737-755, 1997.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the M-by-N matrix to be factored.
          On exit, the factors L and U from the factorization
          A-S=L*U; the unit diagonal elements of L are not stored.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]D
          D is DOUBLE PRECISION array, dimension min(M,N)
          The diagonal elements of the diagonal M-by-N sign matrix S,
          D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can
          be only plus or minus one.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
 November 2019, Igor Kozachenko,
                Computer Science Division,
                University of California, Berkeley

Definition at line 166 of file dlaorhr_col_getrfnp2.f.

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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
recursive subroutine dlaorhr_col_getrfnp2(M, N, A, LDA, D, INFO)
DLAORHR_COL_GETRFNP2
Here is the call graph for this function:
Here is the caller graph for this function: