LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dgetrf.f
Go to the documentation of this file.
1 C> \brief \b DGETRF VARIANT: iterative version of Sivan Toledo's recursive LU algorithm
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER INFO, LDA, M, N
15 * ..
16 * .. Array Arguments ..
17 * INTEGER IPIV( * )
18 * DOUBLE PRECISION A( LDA, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 C>\details \b Purpose:
25 C>\verbatim
26 C>
27 C> DGETRF computes an LU factorization of a general M-by-N matrix A
28 C> using partial pivoting with row interchanges.
29 C>
30 C> The factorization has the form
31 C> A = P * L * U
32 C> where P is a permutation matrix, L is lower triangular with unit
33 C> diagonal elements (lower trapezoidal if m > n), and U is upper
34 C> triangular (upper trapezoidal if m < n).
35 C>
36 C> This code implements an iterative version of Sivan Toledo's recursive
37 C> LU algorithm[1]. For square matrices, this iterative versions should
38 C> be within a factor of two of the optimum number of memory transfers.
39 C>
40 C> The pattern is as follows, with the large blocks of U being updated
41 C> in one call to DTRSM, and the dotted lines denoting sections that
42 C> have had all pending permutations applied:
43 C>
44 C> 1 2 3 4 5 6 7 8
45 C> +-+-+---+-------+------
46 C> | |1| | |
47 C> |.+-+ 2 | |
48 C> | | | | |
49 C> |.|.+-+-+ 4 |
50 C> | | | |1| |
51 C> | | |.+-+ |
52 C> | | | | | |
53 C> |.|.|.|.+-+-+---+ 8
54 C> | | | | | |1| |
55 C> | | | | |.+-+ 2 |
56 C> | | | | | | | |
57 C> | | | | |.|.+-+-+
58 C> | | | | | | | |1|
59 C> | | | | | | |.+-+
60 C> | | | | | | | | |
61 C> |.|.|.|.|.|.|.|.+-----
62 C> | | | | | | | | |
63 C>
64 C> The 1-2-1-4-1-2-1-8-... pattern is the position of the last 1 bit in
65 C> the binary expansion of the current column. Each Schur update is
66 C> applied as soon as the necessary portion of U is available.
67 C>
68 C> [1] Toledo, S. 1997. Locality of Reference in LU Decomposition with
69 C> Partial Pivoting. SIAM J. Matrix Anal. Appl. 18, 4 (Oct. 1997),
70 C> 1065-1081. http://dx.doi.org/10.1137/S0895479896297744
71 C>
72 C>\endverbatim
73 *
74 * Arguments:
75 * ==========
76 *
77 C> \param[in] M
78 C> \verbatim
79 C> M is INTEGER
80 C> The number of rows of the matrix A. M >= 0.
81 C> \endverbatim
82 C>
83 C> \param[in] N
84 C> \verbatim
85 C> N is INTEGER
86 C> The number of columns of the matrix A. N >= 0.
87 C> \endverbatim
88 C>
89 C> \param[in,out] A
90 C> \verbatim
91 C> A is DOUBLE PRECISION array, dimension (LDA,N)
92 C> On entry, the M-by-N matrix to be factored.
93 C> On exit, the factors L and U from the factorization
94 C> A = P*L*U; the unit diagonal elements of L are not stored.
95 C> \endverbatim
96 C>
97 C> \param[in] LDA
98 C> \verbatim
99 C> LDA is INTEGER
100 C> The leading dimension of the array A. LDA >= max(1,M).
101 C> \endverbatim
102 C>
103 C> \param[out] IPIV
104 C> \verbatim
105 C> IPIV is INTEGER array, dimension (min(M,N))
106 C> The pivot indices; for 1 <= i <= min(M,N), row i of the
107 C> matrix was interchanged with row IPIV(i).
108 C> \endverbatim
109 C>
110 C> \param[out] INFO
111 C> \verbatim
112 C> INFO is INTEGER
113 C> = 0: successful exit
114 C> < 0: if INFO = -i, the i-th argument had an illegal value
115 C> > 0: if INFO = i, U(i,i) is exactly zero. The factorization
116 C> has been completed, but the factor U is exactly
117 C> singular, and division by zero will occur if it is used
118 C> to solve a system of equations.
119 C> \endverbatim
120 C>
121 *
122 * Authors:
123 * ========
124 *
125 C> \author Univ. of Tennessee
126 C> \author Univ. of California Berkeley
127 C> \author Univ. of Colorado Denver
128 C> \author NAG Ltd.
129 *
130 C> \date December 2016
131 *
132 C> \ingroup variantsGEcomputational
133 *
134 * =====================================================================
135  SUBROUTINE dgetrf( M, N, A, LDA, IPIV, INFO )
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 *
141 * .. Scalar Arguments ..
142  INTEGER INFO, LDA, M, N
143 * ..
144 * .. Array Arguments ..
145  INTEGER IPIV( * )
146  DOUBLE PRECISION A( LDA, * )
147 * ..
148 *
149 * =====================================================================
150 *
151 * .. Parameters ..
152  DOUBLE PRECISION ONE, ZERO, NEGONE
153  parameter( one = 1.0d+0, zero = 0.0d+0 )
154  parameter( negone = -1.0d+0 )
155 * ..
156 * .. Local Scalars ..
157  DOUBLE PRECISION SFMIN, TMP
158  INTEGER I, J, JP, NSTEP, NTOPIV, NPIVED, KAHEAD
159  INTEGER KSTART, IPIVSTART, JPIVSTART, KCOLS
160 * ..
161 * .. External Functions ..
162  DOUBLE PRECISION DLAMCH
163  INTEGER IDAMAX
164  LOGICAL DISNAN
165  EXTERNAL dlamch, idamax, disnan
166 * ..
167 * .. External Subroutines ..
168  EXTERNAL dtrsm, dscal, xerbla, dlaswp
169 * ..
170 * .. Intrinsic Functions ..
171  INTRINSIC max, min, iand
172 * ..
173 * .. Executable Statements ..
174 *
175 * Test the input parameters.
176 *
177  info = 0
178  IF( m.LT.0 ) THEN
179  info = -1
180  ELSE IF( n.LT.0 ) THEN
181  info = -2
182  ELSE IF( lda.LT.max( 1, m ) ) THEN
183  info = -4
184  END IF
185  IF( info.NE.0 ) THEN
186  CALL xerbla( 'DGETRF', -info )
187  RETURN
188  END IF
189 *
190 * Quick return if possible
191 *
192  IF( m.EQ.0 .OR. n.EQ.0 )
193  $ RETURN
194 *
195 * Compute machine safe minimum
196 *
197  sfmin = dlamch( 'S' )
198 *
199  nstep = min( m, n )
200  DO j = 1, nstep
201  kahead = iand( j, -j )
202  kstart = j + 1 - kahead
203  kcols = min( kahead, m-j )
204 *
205 * Find pivot.
206 *
207  jp = j - 1 + idamax( m-j+1, a( j, j ), 1 )
208  ipiv( j ) = jp
209 
210 * Permute just this column.
211  IF (jp .NE. j) THEN
212  tmp = a( j, j )
213  a( j, j ) = a( jp, j )
214  a( jp, j ) = tmp
215  END IF
216 
217 * Apply pending permutations to L
218  ntopiv = 1
219  ipivstart = j
220  jpivstart = j - ntopiv
221  DO WHILE ( ntopiv .LT. kahead )
222  CALL dlaswp( ntopiv, a( 1, jpivstart ), lda, ipivstart, j,
223  $ ipiv, 1 )
224  ipivstart = ipivstart - ntopiv;
225  ntopiv = ntopiv * 2;
226  jpivstart = jpivstart - ntopiv;
227  END DO
228 
229 * Permute U block to match L
230  CALL dlaswp( kcols, a( 1,j+1 ), lda, kstart, j, ipiv, 1 )
231 
232 * Factor the current column
233  IF( a( j, j ).NE.zero .AND. .NOT.disnan( a( j, j ) ) ) THEN
234  IF( abs(a( j, j )) .GE. sfmin ) THEN
235  CALL dscal( m-j, one / a( j, j ), a( j+1, j ), 1 )
236  ELSE
237  DO i = 1, m-j
238  a( j+i, j ) = a( j+i, j ) / a( j, j )
239  END DO
240  END IF
241  ELSE IF( a( j,j ) .EQ. zero .AND. info .EQ. 0 ) THEN
242  info = j
243  END IF
244 
245 * Solve for U block.
246  CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit', kahead,
247  $ kcols, one, a( kstart, kstart ), lda,
248  $ a( kstart, j+1 ), lda )
249 * Schur complement.
250  CALL dgemm( 'No transpose', 'No transpose', m-j,
251  $ kcols, kahead, negone, a( j+1, kstart ), lda,
252  $ a( kstart, j+1 ), lda, one, a( j+1, j+1 ), lda )
253  END DO
254 
255 * Handle pivot permutations on the way out of the recursion
256  npived = iand( nstep, -nstep )
257  j = nstep - npived
258  DO WHILE ( j .GT. 0 )
259  ntopiv = iand( j, -j )
260  CALL dlaswp( ntopiv, a( 1, j-ntopiv+1 ), lda, j+1, nstep,
261  $ ipiv, 1 )
262  j = j - ntopiv
263  END DO
264 
265 * If short and wide, handle the rest of the columns.
266  IF ( m .LT. n ) THEN
267  CALL dlaswp( n-m, a( 1, m+kcols+1 ), lda, 1, m, ipiv, 1 )
268  CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit', m,
269  $ n-m, one, a, lda, a( 1,m+kcols+1 ), lda )
270  END IF
271 
272  RETURN
273 *
274 * End of DGETRF
275 *
276  END
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
subroutine dgetrf(M, N, A, LDA, IPIV, INFO)
DGETRF
Definition: dgetrf.f:108
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
DLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: dlaswp.f:115