LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zgetrf ( integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)

ZGETRF VARIANT: iterative version of Sivan Toledo's recursive LU algorithm

Purpose:

 ZGETRF computes an LU factorization of a general M-by-N matrix A
 using partial pivoting with row interchanges.

 The factorization has the form
    A = P * L * U
 where P is a permutation matrix, L is lower triangular with unit
 diagonal elements (lower trapezoidal if m > n), and U is upper
 triangular (upper trapezoidal if m < n).

 This code implements an iterative version of Sivan Toledo's recursive
 LU algorithm[1].  For square matrices, this iterative versions should
 be within a factor of two of the optimum number of memory transfers.

 The pattern is as follows, with the large blocks of U being updated
 in one call to DTRSM, and the dotted lines denoting sections that
 have had all pending permutations applied:

  1 2 3 4 5 6 7 8
 +-+-+---+-------+------
 | |1|   |       |
 |.+-+ 2 |       |
 | | |   |       |
 |.|.+-+-+   4   |
 | | | |1|       |
 | | |.+-+       |
 | | | | |       |
 |.|.|.|.+-+-+---+  8
 | | | | | |1|   |
 | | | | |.+-+ 2 |
 | | | | | | |   |
 | | | | |.|.+-+-+
 | | | | | | | |1|
 | | | | | | |.+-+
 | | | | | | | | |
 |.|.|.|.|.|.|.|.+-----
 | | | | | | | | |

 The 1-2-1-4-1-2-1-8-... pattern is the position of the last 1 bit in
 the binary expansion of the current column.  Each Schur update is
 applied as soon as the necessary portion of U is available.

 [1] Toledo, S. 1997. Locality of Reference in LU Decomposition with
 Partial Pivoting. SIAM J. Matrix Anal. Appl. 18, 4 (Oct. 1997),
 1065-1081. http://dx.doi.org/10.1137/S0895479896297744
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 = P*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]IPIV
          IPIV is INTEGER array, dimension (min(M,N))
          The pivot indices; for 1 <= i <= min(M,N), row i of the
          matrix was interchanged with row IPIV(i).
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
                has been completed, but the factor U is exactly
                singular, and division by zero will occur if it is used
                to solve a system of equations.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 136 of file zgetrf.f.

136 *
137 * -- LAPACK computational routine (version 3.X) --
138 * -- LAPACK is a software package provided by Univ. of Tennessee, --
139 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140 * November 2011
141 *
142 * .. Scalar Arguments ..
143  INTEGER info, lda, m, n
144 * ..
145 * .. Array Arguments ..
146  INTEGER ipiv( * )
147  COMPLEX*16 a( lda, * )
148 * ..
149 *
150 * =====================================================================
151 *
152 * .. Parameters ..
153  COMPLEX*16 one, negone
154  DOUBLE PRECISION zero
155  parameter ( one = (1.0d+0, 0.0d+0) )
156  parameter ( negone = (-1.0d+0, 0.0d+0) )
157  parameter ( zero = 0.0d+0 )
158 * ..
159 * .. Local Scalars ..
160  DOUBLE PRECISION sfmin, pivmag
161  COMPLEX*16 tmp
162  INTEGER i, j, jp, nstep, ntopiv, npived, kahead
163  INTEGER kstart, ipivstart, jpivstart, kcols
164 * ..
165 * .. External Functions ..
166  DOUBLE PRECISION dlamch
167  INTEGER izamax
168  LOGICAL disnan
169  EXTERNAL dlamch, izamax, disnan
170 * ..
171 * .. External Subroutines ..
172  EXTERNAL ztrsm, zscal, xerbla, zlaswp
173 * ..
174 * .. Intrinsic Functions ..
175  INTRINSIC max, min, iand, abs
176 * ..
177 * .. Executable Statements ..
178 *
179 * Test the input parameters.
180 *
181  info = 0
182  IF( m.LT.0 ) THEN
183  info = -1
184  ELSE IF( n.LT.0 ) THEN
185  info = -2
186  ELSE IF( lda.LT.max( 1, m ) ) THEN
187  info = -4
188  END IF
189  IF( info.NE.0 ) THEN
190  CALL xerbla( 'ZGETRF', -info )
191  RETURN
192  END IF
193 *
194 * Quick return if possible
195 *
196  IF( m.EQ.0 .OR. n.EQ.0 )
197  $ RETURN
198 *
199 * Compute machine safe minimum
200 *
201  sfmin = dlamch( 'S' )
202 *
203  nstep = min( m, n )
204  DO j = 1, nstep
205  kahead = iand( j, -j )
206  kstart = j + 1 - kahead
207  kcols = min( kahead, m-j )
208 *
209 * Find pivot.
210 *
211  jp = j - 1 + izamax( m-j+1, a( j, j ), 1 )
212  ipiv( j ) = jp
213 
214 ! Permute just this column.
215  IF (jp .NE. j) THEN
216  tmp = a( j, j )
217  a( j, j ) = a( jp, j )
218  a( jp, j ) = tmp
219  END IF
220 
221 ! Apply pending permutations to L
222  ntopiv = 1
223  ipivstart = j
224  jpivstart = j - ntopiv
225  DO WHILE ( ntopiv .LT. kahead )
226  CALL zlaswp( ntopiv, a( 1, jpivstart ), lda, ipivstart, j,
227  $ ipiv, 1 )
228  ipivstart = ipivstart - ntopiv;
229  ntopiv = ntopiv * 2;
230  jpivstart = jpivstart - ntopiv;
231  END DO
232 
233 ! Permute U block to match L
234  CALL zlaswp( kcols, a( 1,j+1 ), lda, kstart, j, ipiv, 1 )
235 
236 ! Factor the current column
237  pivmag = abs( a( j, j ) )
238  IF( pivmag.NE.zero .AND. .NOT.disnan( pivmag ) ) THEN
239  IF( pivmag .GE. sfmin ) THEN
240  CALL zscal( m-j, one / a( j, j ), a( j+1, j ), 1 )
241  ELSE
242  DO i = 1, m-j
243  a( j+i, j ) = a( j+i, j ) / a( j, j )
244  END DO
245  END IF
246  ELSE IF( pivmag .EQ. zero .AND. info .EQ. 0 ) THEN
247  info = j
248  END IF
249 
250 ! Solve for U block.
251  CALL ztrsm( 'Left', 'Lower', 'No transpose', 'Unit', kahead,
252  $ kcols, one, a( kstart, kstart ), lda,
253  $ a( kstart, j+1 ), lda )
254 ! Schur complement.
255  CALL zgemm( 'No transpose', 'No transpose', m-j,
256  $ kcols, kahead, negone, a( j+1, kstart ), lda,
257  $ a( kstart, j+1 ), lda, one, a( j+1, j+1 ), lda )
258  END DO
259 
260 ! Handle pivot permutations on the way out of the recursion
261  npived = iand( nstep, -nstep )
262  j = nstep - npived
263  DO WHILE ( j .GT. 0 )
264  ntopiv = iand( j, -j )
265  CALL zlaswp( ntopiv, a( 1, j-ntopiv+1 ), lda, j+1, nstep,
266  $ ipiv, 1 )
267  j = j - ntopiv
268  END DO
269 
270 ! If short and wide, handle the rest of the columns.
271  IF ( m .LT. n ) THEN
272  CALL zlaswp( n-m, a( 1, m+kcols+1 ), lda, 1, m, ipiv, 1 )
273  CALL ztrsm( 'Left', 'Lower', 'No transpose', 'Unit', m,
274  $ n-m, one, a, lda, a( 1,m+kcols+1 ), lda )
275  END IF
276 
277  RETURN
278 *
279 * End of ZGETRF
280 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zlaswp(N, A, LDA, K1, K2, IPIV, INCX)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: zlaswp.f:116
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54

Here is the call graph for this function: