claunhr_col_getrfnp2()

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

CLAUNHR_COL_GETRFNP2

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

Purpose:

 CLAUNHR_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 CUNHR_COL. In CUNHR_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 = [ --&mdash;|--&mdash; ]        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].

 CLAUNHR_COL_GETRFNP2 is called to factorize a block by the blocked
 routine CLAUNHR_COL_GETRFNP, which uses blocked code calling
is self-sufficient and can be used without CLAUNHR_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 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 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.

Date

November 2019

Contributors:

 November 2019, Igor Kozachenko,
                Computer Science Division,
                University of California, Berkeley

Definition at line 169 of file claunhr_col_getrfnp2.f.

      IMPLICIT NONE
*
*  -- LAPACK computational routine (version 3.9.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2019
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, M, N
*     ..
*     .. Array Arguments ..
      COMPLEX         A( LDA, * ), D( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE
      parameter( one = 1.0e+0 )
      COMPLEX            CONE
      parameter( cone = ( 1.0e+0, 0.0e+0 ) )
*     ..
*     .. Local Scalars ..
      REAL               SFMIN
      INTEGER            I, IINFO, N1, N2
      COMPLEX            Z
*     ..
*     .. External Functions ..
      REAL               SLAMCH
      EXTERNAL           slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           cgemm, cscal, ctrsm, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, real, cmplx, aimag, sign, max, min
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      info = 0
      IF( m.LT.0 ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -4
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CLAUNHR_COL_GETRFNP2', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( min( m, n ).EQ.0 )
     $   RETURN
 
      IF ( m.EQ.1 ) THEN
*
*        One row case, (also recursion termination case),
*        use unblocked code
*
*        Transfer the sign
*
         d( 1 ) = cmplx( -sign( one, real( a( 1, 1 ) ) ) )
*
*        Construct the row of U
*
         a( 1, 1 ) = a( 1, 1 ) - d( 1 )
*
      ELSE IF( n.EQ.1 ) THEN
*
*        One column case, (also recursion termination case),
*        use unblocked code
*
*        Transfer the sign
*
         d( 1 ) = cmplx( -sign( one, real( a( 1, 1 ) ) ) )
*
*        Construct the row of U
*
         a( 1, 1 ) = a( 1, 1 ) - d( 1 )
*
*        Scale the elements 2:M of the column
*
*        Determine machine safe minimum
*
         sfmin = slamch('S')
*
*        Construct the subdiagonal elements of L
*
         IF( cabs1( a( 1, 1 ) ) .GE. sfmin ) THEN
            CALL cscal( m-1, cone / a( 1, 1 ), a( 2, 1 ), 1 )
         ELSE
            DO i = 2, m
               a( i, 1 ) = a( i, 1 ) / a( 1, 1 )
            END DO
         END IF
*
      ELSE
*
*        Divide the matrix B into four submatrices
*
         n1 = min( m, n ) / 2
         n2 = n-n1
 
*
*        Factor B11, recursive call
*
         CALL claunhr_col_getrfnp2( n1, n1, a, lda, d, iinfo )
*
*        Solve for B21
*
         CALL ctrsm( 'R', 'U', 'N', 'N', m-n1, n1, cone, a, lda,
     $               a( n1+1, 1 ), lda )
*
*        Solve for B12
*
         CALL ctrsm( 'L', 'L', 'N', 'U', n1, n2, cone, a, lda,
     $               a( 1, n1+1 ), lda )
*
*        Update B22, i.e. compute the Schur complement
*        B22 := B22 - B21*B12
*
         CALL cgemm( 'N', 'N', m-n1, n2, n1, -cone, a( n1+1, 1 ), lda,
     $               a( 1, n1+1 ), lda, cone, a( n1+1, n1+1 ), lda )
*
*        Factor B22, recursive call
*
         CALL claunhr_col_getrfnp2( m-n1, n2, a( n1+1, n1+1 ), lda,
     $                              d( n1+1 ), iinfo )
*
      END IF
      RETURN
*
*     End of CLAUNHR_COL_GETRFNP2
*

Here is the call graph for this function:

Here is the caller graph for this function: