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

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

Purpose:

``` SGETRF 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 STRSM, 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 REAL 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.```
Date
November 2011

Definition at line 136 of file sgetrf.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  REAL a( lda, * )
148 * ..
149 *
150 * =====================================================================
151 *
152 * .. Parameters ..
153  REAL one, zero, negone
154  parameter ( one = 1.0e+0, zero = 0.0e+0 )
155  parameter ( negone = -1.0e+0 )
156 * ..
157 * .. Local Scalars ..
158  REAL sfmin, tmp
159  INTEGER i, j, jp, nstep, ntopiv, npived, kahead
160  INTEGER kstart, ipivstart, jpivstart, kcols
161 * ..
162 * .. External Functions ..
163  REAL slamch
164  INTEGER isamax
165  LOGICAL sisnan
166  EXTERNAL slamch, isamax, sisnan
167 * ..
168 * .. External Subroutines ..
169  EXTERNAL strsm, sscal, xerbla, slaswp
170 * ..
171 * .. Intrinsic Functions ..
172  INTRINSIC max, min, iand
173 * ..
174 * .. Executable Statements ..
175 *
176 * Test the input parameters.
177 *
178  info = 0
179  IF( m.LT.0 ) THEN
180  info = -1
181  ELSE IF( n.LT.0 ) THEN
182  info = -2
183  ELSE IF( lda.LT.max( 1, m ) ) THEN
184  info = -4
185  END IF
186  IF( info.NE.0 ) THEN
187  CALL xerbla( 'SGETRF', -info )
188  RETURN
189  END IF
190 *
191 * Quick return if possible
192 *
193  IF( m.EQ.0 .OR. n.EQ.0 )
194  \$ RETURN
195 *
196 * Compute machine safe minimum
197 *
198  sfmin = slamch( 'S' )
199 *
200  nstep = min( m, n )
201  DO j = 1, nstep
202  kahead = iand( j, -j )
203  kstart = j + 1 - kahead
204  kcols = min( kahead, m-j )
205 *
206 * Find pivot.
207 *
208  jp = j - 1 + isamax( m-j+1, a( j, j ), 1 )
209  ipiv( j ) = jp
210
211 ! Permute just this column.
212  IF (jp .NE. j) THEN
213  tmp = a( j, j )
214  a( j, j ) = a( jp, j )
215  a( jp, j ) = tmp
216  END IF
217
218 ! Apply pending permutations to L
219  ntopiv = 1
220  ipivstart = j
221  jpivstart = j - ntopiv
222  DO WHILE ( ntopiv .LT. kahead )
223  CALL slaswp( ntopiv, a( 1, jpivstart ), lda, ipivstart, j,
224  \$ ipiv, 1 )
225  ipivstart = ipivstart - ntopiv;
226  ntopiv = ntopiv * 2;
227  jpivstart = jpivstart - ntopiv;
228  END DO
229
230 ! Permute U block to match L
231  CALL slaswp( kcols, a( 1,j+1 ), lda, kstart, j, ipiv, 1 )
232
233 ! Factor the current column
234  IF( a( j, j ).NE.zero .AND. .NOT.sisnan( a( j, j ) ) ) THEN
235  IF( abs(a( j, j )) .GE. sfmin ) THEN
236  CALL sscal( m-j, one / a( j, j ), a( j+1, j ), 1 )
237  ELSE
238  DO i = 1, m-j
239  a( j+i, j ) = a( j+i, j ) / a( j, j )
240  END DO
241  END IF
242  ELSE IF( a( j,j ) .EQ. zero .AND. info .EQ. 0 ) THEN
243  info = j
244  END IF
245
246 ! Solve for U block.
247  CALL strsm( 'Left', 'Lower', 'No transpose', 'Unit', kahead,
248  \$ kcols, one, a( kstart, kstart ), lda,
249  \$ a( kstart, j+1 ), lda )
250 ! Schur complement.
251  CALL sgemm( 'No transpose', 'No transpose', m-j,
252  \$ kcols, kahead, negone, a( j+1, kstart ), lda,
253  \$ a( kstart, j+1 ), lda, one, a( j+1, j+1 ), lda )
254  END DO
255
256 ! Handle pivot permutations on the way out of the recursion
257  npived = iand( nstep, -nstep )
258  j = nstep - npived
259  DO WHILE ( j .GT. 0 )
260  ntopiv = iand( j, -j )
261  CALL slaswp( ntopiv, a( 1, j-ntopiv+1 ), lda, j+1, nstep,
262  \$ ipiv, 1 )
263  j = j - ntopiv
264  END DO
265
266 ! If short and wide, handle the rest of the columns.
267  IF ( m .LT. n ) THEN
268  CALL slaswp( n-m, a( 1, m+kcols+1 ), lda, 1, m, ipiv, 1 )
269  CALL strsm( 'Left', 'Lower', 'No transpose', 'Unit', m,
270  \$ n-m, one, a, lda, a( 1,m+kcols+1 ), lda )
271  END IF
272
273  RETURN
274 *
275 * End of SGETRF
276 *
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
logical function sisnan(SIN)
SISNAN tests input for NaN.
Definition: sisnan.f:61
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slaswp(N, A, LDA, K1, K2, IPIV, INCX)
SLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: slaswp.f:116
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69

Here is the call graph for this function: