LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zlaunhr_col_getrfnp2()

recursive subroutine zlaunhr_col_getrfnp2 ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  D,
integer  INFO 
)

ZLAUNHR_COL_GETRFNP2

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

Purpose:
 ZLAUNHR_COL_GETRFNP2 computes the modified LU factorization without
 pivoting of a complex 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 ZUNHR_COL. In ZUNHR_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].

 ZLAUNHR_COL_GETRFNP2 is called to factorize a block by the blocked
 routine ZLAUNHR_COL_GETRFNP, which uses blocked code calling
 Level 3 BLAS to update the submatrix. However, ZLAUNHR_COL_GETRFNP2
 is self-sufficient and can be used without ZLAUNHR_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 COMPLEX*16 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 COMPLEX*16 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 ( +1.0, 0.0 ) or (-1.0, 0.0 ).
[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 zlaunhr_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  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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:180
recursive subroutine zlaunhr_col_getrfnp2(M, N, A, LDA, D, INFO)
ZLAUNHR_COL_GETRFNP2
Here is the call graph for this function:
Here is the caller graph for this function: