OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
Variants Computational routines
LAPACK » General Matrices

Functions

subroutine cgetrf (m, n, a, lda, ipiv, info)
 CGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
subroutine dgetrf (m, n, a, lda, ipiv, info)
 DGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
subroutine sgetrf (m, n, a, lda, ipiv, info)
 SGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
subroutine zgetrf (m, n, a, lda, ipiv, info)
 ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
subroutine cgeqrf (m, n, a, lda, tau, work, lwork, info)
 CGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.
subroutine dgeqrf (m, n, a, lda, tau, work, lwork, info)
 DGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.
subroutine sgeqrf (m, n, a, lda, tau, work, lwork, info)
 SGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.
subroutine zgeqrf (m, n, a, lda, tau, work, lwork, info)
 ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.

Detailed Description

This is the group of Variants Computational routines

Function Documentation

◆ cgeqrf()

subroutine cgeqrf ( integer m,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( * ) tau,
complex, dimension( * ) work,
integer lwork,
integer info )

CGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.

Purpose: 

!>
!> CGEQRF computes a QR factorization of a real M-by-N matrix A:
!> A = Q * R.
!>
!> This is the left-looking Level 3 BLAS version of the algorithm.
!>
!>
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 A.
!>          On exit, the elements on and above the diagonal of the array
!>          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
!>          upper triangular if m >= n); the elements below the diagonal,
!>          with the array TAU, represent the orthogonal matrix Q as a
!>          product of min(m,n) elementary reflectors (see Further
!>          Details).
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[out]TAU
!>          TAU is COMPLEX array, dimension (min(M,N))
!>          The scalar factors of the elementary reflectors (see Further
!>          Details).
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!> 
!>          The dimension of the array WORK. The dimension can be divided into three parts.
!> 
!>          1) The part for the triangular factor T. If the very last T is not bigger
!>             than any of the rest, then this part is NB x ceiling(K/NB), otherwise,
!>             NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T
!> 
!>          2) The part for the very last T when T is bigger than any of the rest T.
!>             The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB,
!>             where K = min(M,N), NX is calculated by
!>                   NX = MAX( 0, ILAENV( 3, 'CGEQRF', ' ', M, N, -1, -1 ) )
!> 
!>          3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB)
!> 
!>          So LWORK = part1 + part2 + part3
!> 
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[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
December 2016

Further Details 

!>
!>  The matrix Q is represented as a product of elementary reflectors
!>
!>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
!>
!>  Each H(i) has the form
!>
!>     H(i) = I - tau * v * v'
!>
!>  where tau is a real scalar, and v is a real vector with
!>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
!>  and tau in TAU(i).
!>
!>

Definition at line 150 of file cgeqrf.f.

151*
152* -- LAPACK computational routine --
153* -- LAPACK is a software package provided by Univ. of Tennessee, --
154* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155*
156* .. Scalar Arguments ..
157 INTEGER INFO, LDA, LWORK, M, N
158* ..
159* .. Array Arguments ..
160 COMPLEX A( LDA, * ), TAU( * ), WORK( * )
161* ..
162*
163* =====================================================================
164*
165* .. Local Scalars ..
166 LOGICAL LQUERY
167 INTEGER I, IB, IINFO, IWS, J, K, LWKOPT, NB,
168 $ NBMIN, NX, LBWORK, NT, LLWORK
169* ..
170* .. External Subroutines ..
171 EXTERNAL cgeqr2, clarfb, clarft, xerbla
172* ..
173* .. Intrinsic Functions ..
174 INTRINSIC max, min
175* ..
176* .. External Functions ..
177 INTEGER ILAENV
178 REAL SCEIL
179 EXTERNAL ilaenv, sceil
180* ..
181* .. Executable Statements ..
182
183 info = 0
184 nbmin = 2
185 nx = 0
186 iws = n
187 k = min( m, n )
188 nb = ilaenv( 1, 'CGEQRF', ' ', m, n, -1, -1 )
189
190 IF( nb.GT.1 .AND. nb.LT.k ) THEN
191*
192* Determine when to cross over from blocked to unblocked code.
193*
194 nx = max( 0, ilaenv( 3, 'CGEQRF', ' ', m, n, -1, -1 ) )
195 END IF
196*
197* Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.:
198*
199* NB=3 2NB=6 K=10
200* | | |
201* 1--2--3--4--5--6--7--8--9--10
202* | \________/
203* K-NX=5 NT=4
204*
205* So here 4 x 4 is the last T stored in the workspace
206*
207 nt = k-sceil(real(k-nx)/real(nb))*nb
208
209*
210* optimal workspace = space for dlarfb + space for normal T's + space for the last T
211*
212 llwork = max(max((n-m)*k, (n-m)*nb), max(k*nb, nb*nb))
213 llwork = sceil(real(llwork)/real(nb))
214
215 IF ( nt.GT.nb ) THEN
216
217 lbwork = k-nt
218*
219* Optimal workspace for dlarfb = MAX(1,N)*NT
220*
221 lwkopt = (lbwork+llwork)*nb
222 work( 1 ) = (lwkopt+nt*nt)
223
224 ELSE
225
226 lbwork = sceil(real(k)/real(nb))*nb
227 lwkopt = (lbwork+llwork-nb)*nb
228 work( 1 ) = lwkopt
229
230 END IF
231
232*
233* Test the input arguments
234*
235 lquery = ( lwork.EQ.-1 )
236 IF( m.LT.0 ) THEN
237 info = -1
238 ELSE IF( n.LT.0 ) THEN
239 info = -2
240 ELSE IF( lda.LT.max( 1, m ) ) THEN
241 info = -4
242 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
243 info = -7
244 END IF
245 IF( info.NE.0 ) THEN
246 CALL xerbla( 'CGEQRF', -info )
247 RETURN
248 ELSE IF( lquery ) THEN
249 RETURN
250 END IF
251*
252* Quick return if possible
253*
254 IF( k.EQ.0 ) THEN
255 work( 1 ) = 1
256 RETURN
257 END IF
258*
259 IF( nb.GT.1 .AND. nb.LT.k ) THEN
260
261 IF( nx.LT.k ) THEN
262*
263* Determine if workspace is large enough for blocked code.
264*
265 IF ( nt.LE.nb ) THEN
266 iws = (lbwork+llwork-nb)*nb
267 ELSE
268 iws = (lbwork+llwork)*nb+nt*nt
269 END IF
270
271 IF( lwork.LT.iws ) THEN
272*
273* Not enough workspace to use optimal NB: reduce NB and
274* determine the minimum value of NB.
275*
276 IF ( nt.LE.nb ) THEN
277 nb = lwork / (llwork+(lbwork-nb))
278 ELSE
279 nb = (lwork-nt*nt)/(lbwork+llwork)
280 END IF
281
282 nbmin = max( 2, ilaenv( 2, 'CGEQRF', ' ', m, n, -1,
283 $ -1 ) )
284 END IF
285 END IF
286 END IF
287*
288 IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
289*
290* Use blocked code initially
291*
292 DO 10 i = 1, k - nx, nb
293 ib = min( k-i+1, nb )
294*
295* Update the current column using old T's
296*
297 DO 20 j = 1, i - nb, nb
298*
299* Apply H' to A(J:M,I:I+IB-1) from the left
300*
301 CALL clarfb( 'Left', 'Transpose', 'Forward',
302 $ 'Columnwise', m-j+1, ib, nb,
303 $ a( j, j ), lda, work(j), lbwork,
304 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
305 $ ib)
306
30720 CONTINUE
308*
309* Compute the QR factorization of the current block
310* A(I:M,I:I+IB-1)
311*
312 CALL cgeqr2( m-i+1, ib, a( i, i ), lda, tau( i ),
313 $ work(lbwork*nb+nt*nt+1), iinfo )
314
315 IF( i+ib.LE.n ) THEN
316*
317* Form the triangular factor of the block reflector
318* H = H(i) H(i+1) . . . H(i+ib-1)
319*
320 CALL clarft( 'Forward', 'Columnwise', m-i+1, ib,
321 $ a( i, i ), lda, tau( i ),
322 $ work(i), lbwork )
323*
324 END IF
325 10 CONTINUE
326 ELSE
327 i = 1
328 END IF
329*
330* Use unblocked code to factor the last or only block.
331*
332 IF( i.LE.k ) THEN
333
334 IF ( i .NE. 1 ) THEN
335
336 DO 30 j = 1, i - nb, nb
337*
338* Apply H' to A(J:M,I:K) from the left
339*
340 CALL clarfb( 'Left', 'Transpose', 'Forward',
341 $ 'Columnwise', m-j+1, k-i+1, nb,
342 $ a( j, j ), lda, work(j), lbwork,
343 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
344 $ k-i+1)
34530 CONTINUE
346
347 CALL cgeqr2( m-i+1, k-i+1, a( i, i ), lda, tau( i ),
348 $ work(lbwork*nb+nt*nt+1),iinfo )
349
350 ELSE
351*
352* Use unblocked code to factor the last or only block.
353*
354 CALL cgeqr2( m-i+1, n-i+1, a( i, i ), lda, tau( i ),
355 $ work,iinfo )
356
357 END IF
358 END IF
359
360
361*
362* Apply update to the column M+1:N when N > M
363*
364 IF ( m.LT.n .AND. i.NE.1) THEN
365*
366* Form the last triangular factor of the block reflector
367* H = H(i) H(i+1) . . . H(i+ib-1)
368*
369 IF ( nt .LE. nb ) THEN
370 CALL clarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
371 $ a( i, i ), lda, tau( i ), work(i), lbwork )
372 ELSE
373 CALL clarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
374 $ a( i, i ), lda, tau( i ),
375 $ work(lbwork*nb+1), nt )
376 END IF
377
378*
379* Apply H' to A(1:M,M+1:N) from the left
380*
381 DO 40 j = 1, k-nx, nb
382
383 ib = min( k-j+1, nb )
384
385 CALL clarfb( 'Left', 'Transpose', 'Forward',
386 $ 'Columnwise', m-j+1, n-m, ib,
387 $ a( j, j ), lda, work(j), lbwork,
388 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
389 $ n-m)
390
39140 CONTINUE
392
393 IF ( nt.LE.nb ) THEN
394 CALL clarfb( 'Left', 'Transpose', 'Forward',
395 $ 'Columnwise', m-j+1, n-m, k-j+1,
396 $ a( j, j ), lda, work(j), lbwork,
397 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
398 $ n-m)
399 ELSE
400 CALL clarfb( 'Left', 'Transpose', 'Forward',
401 $ 'Columnwise', m-j+1, n-m, k-j+1,
402 $ a( j, j ), lda,
403 $ work(lbwork*nb+1),
404 $ nt, a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
405 $ n-m)
406 END IF
407
408 END IF
409
410 work( 1 ) = iws
411 RETURN
412*
413* End of CGEQRF
414*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
cgeqr2
subroutine cgeqr2(m, n, a, lda, tau, work, info)
CGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition cgeqr2.f:130
clarfb
subroutine clarfb(side, trans, direct, storev, m, n, k, v, ldv, t, ldt, c, ldc, work, ldwork)
CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
Definition clarfb.f:197
clarft
subroutine clarft(direct, storev, n, k, v, ldv, tau, t, ldt)
CLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition clarft.f:163
sceil
real function sceil(a)
SCEIL
Definition sceil.f:60
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ cgetrf()

subroutine cgetrf ( integer m,
integer n,
complex, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer info )

CGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.

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

CGETRF VARIANT: left-looking Level 3 BLAS version of the algorithm.

Purpose: 

!>
!> CGETRF 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 is the Crout Level 3 BLAS version of the algorithm.
!>
!>
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 = 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
December 2016

Purpose: 

!>
!> CGETRF 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 is the left-looking Level 3 BLAS version of the algorithm.
!>
!>
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 = 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
December 2016

Purpose: 

!>
!> CGETRF 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 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
December 2016

Definition at line 101 of file cgetrf.f.

102*
103* -- LAPACK computational routine --
104* -- LAPACK is a software package provided by Univ. of Tennessee, --
105* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
106*
107* .. Scalar Arguments ..
108 INTEGER INFO, LDA, M, N
109* ..
110* .. Array Arguments ..
111 INTEGER IPIV( * )
112 COMPLEX A( LDA, * )
113* ..
114*
115* =====================================================================
116*
117* .. Parameters ..
118 COMPLEX ONE
119 parameter( one = ( 1.0e+0, 0.0e+0 ) )
120* ..
121* .. Local Scalars ..
122 INTEGER I, IINFO, J, JB, NB
123* ..
124* .. External Subroutines ..
125 EXTERNAL cgemm, cgetf2, claswp, ctrsm, xerbla
126* ..
127* .. External Functions ..
128 INTEGER ILAENV
129 EXTERNAL ilaenv
130* ..
131* .. Intrinsic Functions ..
132 INTRINSIC max, min
133* ..
134* .. Executable Statements ..
135*
136* Test the input parameters.
137*
138 info = 0
139 IF( m.LT.0 ) THEN
140 info = -1
141 ELSE IF( n.LT.0 ) THEN
142 info = -2
143 ELSE IF( lda.LT.max( 1, m ) ) THEN
144 info = -4
145 END IF
146 IF( info.NE.0 ) THEN
147 CALL xerbla( 'CGETRF', -info )
148 RETURN
149 END IF
150*
151* Quick return if possible
152*
153 IF( m.EQ.0 .OR. n.EQ.0 )
154 $ RETURN
155*
156* Determine the block size for this environment.
157*
158 nb = ilaenv( 1, 'CGETRF', ' ', m, n, -1, -1 )
159 IF( nb.LE.1 .OR. nb.GE.min( m, n ) ) THEN
160*
161* Use unblocked code.
162*
163 CALL cgetf2( m, n, a, lda, ipiv, info )
164 ELSE
165*
166* Use blocked code.
167*
168 DO 20 j = 1, min( m, n ), nb
169 jb = min( min( m, n )-j+1, nb )
170*
171* Update current block.
172*
173 CALL cgemm( 'No transpose', 'No transpose',
174 $ m-j+1, jb, j-1, -one,
175 $ a( j, 1 ), lda, a( 1, j ), lda, one,
176 $ a( j, j ), lda )
177
178*
179* Factor diagonal and subdiagonal blocks and test for exact
180* singularity.
181*
182 CALL cgetf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo )
183*
184* Adjust INFO and the pivot indices.
185*
186 IF( info.EQ.0 .AND. iinfo.GT.0 )
187 $ info = iinfo + j - 1
188 DO 10 i = j, min( m, j+jb-1 )
189 ipiv( i ) = j - 1 + ipiv( i )
190 10 CONTINUE
191*
192* Apply interchanges to column 1:J-1
193*
194 CALL claswp( j-1, a, lda, j, j+jb-1, ipiv, 1 )
195*
196 IF ( j+jb.LE.n ) THEN
197*
198* Apply interchanges to column J+JB:N
199*
200 CALL claswp( n-j-jb+1, a( 1, j+jb ), lda, j, j+jb-1,
201 $ ipiv, 1 )
202*
203 CALL cgemm( 'No transpose', 'No transpose',
204 $ jb, n-j-jb+1, j-1, -one,
205 $ a( j, 1 ), lda, a( 1, j+jb ), lda, one,
206 $ a( j, j+jb ), lda )
207*
208* Compute block row of U.
209*
210 CALL ctrsm( 'Left', 'Lower', 'No transpose', 'Unit',
211 $ jb, n-j-jb+1, one, a( j, j ), lda,
212 $ a( j, j+jb ), lda )
213 END IF
214
215 20 CONTINUE
216
217 END IF
218 RETURN
219*
220* End of CGETRF
221*
cgetf2
subroutine cgetf2(m, n, a, lda, ipiv, info)
CGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row inter...
Definition cgetf2.f:108
claswp
subroutine claswp(n, a, lda, k1, k2, ipiv, incx)
CLASWP performs a series of row interchanges on a general rectangular matrix.
Definition claswp.f:115
ctrsm
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187

◆ dgeqrf()

subroutine dgeqrf ( integer m,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) tau,
double precision, dimension( * ) work,
integer lwork,
integer info )

DGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.

Purpose: 

!>
!> DGEQRF computes a QR factorization of a real M-by-N matrix A:
!> A = Q * R.
!>
!> This is the left-looking Level 3 BLAS version of the algorithm.
!>
!>
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 DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit, the elements on and above the diagonal of the array
!>          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
!>          upper triangular if m >= n); the elements below the diagonal,
!>          with the array TAU, represent the orthogonal matrix Q as a
!>          product of min(m,n) elementary reflectors (see Further
!>          Details).
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[out]TAU
!>          TAU is DOUBLE PRECISION array, dimension (min(M,N))
!>          The scalar factors of the elementary reflectors (see Further
!>          Details).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!> 
!>          The dimension of the array WORK. The dimension can be divided into three parts.
!> 
!>          1) The part for the triangular factor T. If the very last T is not bigger
!>             than any of the rest, then this part is NB x ceiling(K/NB), otherwise,
!>             NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T
!> 
!>          2) The part for the very last T when T is bigger than any of the rest T.
!>             The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB,
!>             where K = min(M,N), NX is calculated by
!>                   NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) )
!> 
!>          3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB)
!> 
!>          So LWORK = part1 + part2 + part3
!> 
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[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
December 2016

Further Details 

!>
!>  The matrix Q is represented as a product of elementary reflectors
!>
!>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
!>
!>  Each H(i) has the form
!>
!>     H(i) = I - tau * v * v'
!>
!>  where tau is a real scalar, and v is a real vector with
!>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
!>  and tau in TAU(i).
!>
!>

Definition at line 150 of file dgeqrf.f.

151*
152* -- LAPACK computational routine --
153* -- LAPACK is a software package provided by Univ. of Tennessee, --
154* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155*
156* .. Scalar Arguments ..
157 INTEGER INFO, LDA, LWORK, M, N
158* ..
159* .. Array Arguments ..
160 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
161* ..
162*
163* =====================================================================
164*
165* .. Local Scalars ..
166 LOGICAL LQUERY
167 INTEGER I, IB, IINFO, IWS, J, K, LWKOPT, NB,
168 $ NBMIN, NX, LBWORK, NT, LLWORK
169* ..
170* .. External Subroutines ..
171 EXTERNAL dgeqr2, dlarfb, dlarft, xerbla
172* ..
173* .. Intrinsic Functions ..
174 INTRINSIC max, min
175* ..
176* .. External Functions ..
177 INTEGER ILAENV
178 REAL SCEIL
179 EXTERNAL ilaenv, sceil
180* ..
181* .. Executable Statements ..
182
183 info = 0
184 nbmin = 2
185 nx = 0
186 iws = n
187 k = min( m, n )
188 nb = ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
189
190 IF( nb.GT.1 .AND. nb.LT.k ) THEN
191*
192* Determine when to cross over from blocked to unblocked code.
193*
194 nx = max( 0, ilaenv( 3, 'DGEQRF', ' ', m, n, -1, -1 ) )
195 END IF
196*
197* Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.:
198*
199* NB=3 2NB=6 K=10
200* | | |
201* 1--2--3--4--5--6--7--8--9--10
202* | \________/
203* K-NX=5 NT=4
204*
205* So here 4 x 4 is the last T stored in the workspace
206*
207 nt = k-sceil(real(k-nx)/real(nb))*nb
208
209*
210* optimal workspace = space for dlarfb + space for normal T's + space for the last T
211*
212 llwork = max(max((n-m)*k, (n-m)*nb), max(k*nb, nb*nb))
213 llwork = sceil(real(llwork)/real(nb))
214
215 IF ( nt.GT.nb ) THEN
216
217 lbwork = k-nt
218*
219* Optimal workspace for dlarfb = MAX(1,N)*NT
220*
221 lwkopt = (lbwork+llwork)*nb
222 work( 1 ) = (lwkopt+nt*nt)
223
224 ELSE
225
226 lbwork = sceil(real(k)/real(nb))*nb
227 lwkopt = (lbwork+llwork-nb)*nb
228 work( 1 ) = lwkopt
229
230 END IF
231
232*
233* Test the input arguments
234*
235 lquery = ( lwork.EQ.-1 )
236 IF( m.LT.0 ) THEN
237 info = -1
238 ELSE IF( n.LT.0 ) THEN
239 info = -2
240 ELSE IF( lda.LT.max( 1, m ) ) THEN
241 info = -4
242 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
243 info = -7
244 END IF
245 IF( info.NE.0 ) THEN
246 CALL xerbla( 'DGEQRF', -info )
247 RETURN
248 ELSE IF( lquery ) THEN
249 RETURN
250 END IF
251*
252* Quick return if possible
253*
254 IF( k.EQ.0 ) THEN
255 work( 1 ) = 1
256 RETURN
257 END IF
258*
259 IF( nb.GT.1 .AND. nb.LT.k ) THEN
260
261 IF( nx.LT.k ) THEN
262*
263* Determine if workspace is large enough for blocked code.
264*
265 IF ( nt.LE.nb ) THEN
266 iws = (lbwork+llwork-nb)*nb
267 ELSE
268 iws = (lbwork+llwork)*nb+nt*nt
269 END IF
270
271 IF( lwork.LT.iws ) THEN
272*
273* Not enough workspace to use optimal NB: reduce NB and
274* determine the minimum value of NB.
275*
276 IF ( nt.LE.nb ) THEN
277 nb = lwork / (llwork+(lbwork-nb))
278 ELSE
279 nb = (lwork-nt*nt)/(lbwork+llwork)
280 END IF
281
282 nbmin = max( 2, ilaenv( 2, 'DGEQRF', ' ', m, n, -1,
283 $ -1 ) )
284 END IF
285 END IF
286 END IF
287*
288 IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
289*
290* Use blocked code initially
291*
292 DO 10 i = 1, k - nx, nb
293 ib = min( k-i+1, nb )
294*
295* Update the current column using old T's
296*
297 DO 20 j = 1, i - nb, nb
298*
299* Apply H' to A(J:M,I:I+IB-1) from the left
300*
301 CALL dlarfb( 'Left', 'Transpose', 'Forward',
302 $ 'Columnwise', m-j+1, ib, nb,
303 $ a( j, j ), lda, work(j), lbwork,
304 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
305 $ ib)
306
30720 CONTINUE
308*
309* Compute the QR factorization of the current block
310* A(I:M,I:I+IB-1)
311*
312 CALL dgeqr2( m-i+1, ib, a( i, i ), lda, tau( i ),
313 $ work(lbwork*nb+nt*nt+1), iinfo )
314
315 IF( i+ib.LE.n ) THEN
316*
317* Form the triangular factor of the block reflector
318* H = H(i) H(i+1) . . . H(i+ib-1)
319*
320 CALL dlarft( 'Forward', 'Columnwise', m-i+1, ib,
321 $ a( i, i ), lda, tau( i ),
322 $ work(i), lbwork )
323*
324 END IF
325 10 CONTINUE
326 ELSE
327 i = 1
328 END IF
329*
330* Use unblocked code to factor the last or only block.
331*
332 IF( i.LE.k ) THEN
333
334 IF ( i .NE. 1 ) THEN
335
336 DO 30 j = 1, i - nb, nb
337*
338* Apply H' to A(J:M,I:K) from the left
339*
340 CALL dlarfb( 'Left', 'Transpose', 'Forward',
341 $ 'Columnwise', m-j+1, k-i+1, nb,
342 $ a( j, j ), lda, work(j), lbwork,
343 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
344 $ k-i+1)
34530 CONTINUE
346
347 CALL dgeqr2( m-i+1, k-i+1, a( i, i ), lda, tau( i ),
348 $ work(lbwork*nb+nt*nt+1),iinfo )
349
350 ELSE
351*
352* Use unblocked code to factor the last or only block.
353*
354 CALL dgeqr2( m-i+1, n-i+1, a( i, i ), lda, tau( i ),
355 $ work,iinfo )
356
357 END IF
358 END IF
359
360
361*
362* Apply update to the column M+1:N when N > M
363*
364 IF ( m.LT.n .AND. i.NE.1) THEN
365*
366* Form the last triangular factor of the block reflector
367* H = H(i) H(i+1) . . . H(i+ib-1)
368*
369 IF ( nt .LE. nb ) THEN
370 CALL dlarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
371 $ a( i, i ), lda, tau( i ), work(i), lbwork )
372 ELSE
373 CALL dlarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
374 $ a( i, i ), lda, tau( i ),
375 $ work(lbwork*nb+1), nt )
376 END IF
377
378*
379* Apply H' to A(1:M,M+1:N) from the left
380*
381 DO 40 j = 1, k-nx, nb
382
383 ib = min( k-j+1, nb )
384
385 CALL dlarfb( 'Left', 'Transpose', 'Forward',
386 $ 'Columnwise', m-j+1, n-m, ib,
387 $ a( j, j ), lda, work(j), lbwork,
388 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
389 $ n-m)
390
39140 CONTINUE
392
393 IF ( nt.LE.nb ) THEN
394 CALL dlarfb( 'Left', 'Transpose', 'Forward',
395 $ 'Columnwise', m-j+1, n-m, k-j+1,
396 $ a( j, j ), lda, work(j), lbwork,
397 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
398 $ n-m)
399 ELSE
400 CALL dlarfb( 'Left', 'Transpose', 'Forward',
401 $ 'Columnwise', m-j+1, n-m, k-j+1,
402 $ a( j, j ), lda,
403 $ work(lbwork*nb+1),
404 $ nt, a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
405 $ n-m)
406 END IF
407
408 END IF
409
410 work( 1 ) = iws
411 RETURN
412*
413* End of DGEQRF
414*
dgeqr2
subroutine dgeqr2(m, n, a, lda, tau, work, info)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition dgeqr2.f:130
dlarft
subroutine dlarft(direct, storev, n, k, v, ldv, tau, t, ldt)
DLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition dlarft.f:163
dlarfb
subroutine dlarfb(side, trans, direct, storev, m, n, k, v, ldv, t, ldt, c, ldc, work, ldwork)
DLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition dlarfb.f:197

◆ dgetrf()

subroutine dgetrf ( integer m,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer info )

DGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.

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

DGETRF VARIANT: left-looking Level 3 BLAS version of the algorithm.

Purpose: 

!>
!> DGETRF 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 is the Crout Level 3 BLAS version of the algorithm.
!>
!>
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 DOUBLE PRECISION 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
December 2016

Purpose: 

!>
!> DGETRF 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 is the left-looking Level 3 BLAS version of the algorithm.
!>
!>
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 DOUBLE PRECISION 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
December 2016

Purpose: 

!>
!> DGETRF 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 DOUBLE PRECISION 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
December 2016

Definition at line 101 of file dgetrf.f.

102*
103* -- LAPACK computational routine --
104* -- LAPACK is a software package provided by Univ. of Tennessee, --
105* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
106*
107* .. Scalar Arguments ..
108 INTEGER INFO, LDA, M, N
109* ..
110* .. Array Arguments ..
111 INTEGER IPIV( * )
112 DOUBLE PRECISION A( LDA, * )
113* ..
114*
115* =====================================================================
116*
117* .. Parameters ..
118 DOUBLE PRECISION ONE
119 parameter( one = 1.0d+0 )
120* ..
121* .. Local Scalars ..
122 INTEGER I, IINFO, J, JB, NB
123* ..
124* .. External Subroutines ..
125 EXTERNAL dgemm, dgetf2, dlaswp, dtrsm, xerbla
126* ..
127* .. External Functions ..
128 INTEGER ILAENV
129 EXTERNAL ilaenv
130* ..
131* .. Intrinsic Functions ..
132 INTRINSIC max, min
133* ..
134* .. Executable Statements ..
135*
136* Test the input parameters.
137*
138 info = 0
139 IF( m.LT.0 ) THEN
140 info = -1
141 ELSE IF( n.LT.0 ) THEN
142 info = -2
143 ELSE IF( lda.LT.max( 1, m ) ) THEN
144 info = -4
145 END IF
146 IF( info.NE.0 ) THEN
147 CALL xerbla( 'DGETRF', -info )
148 RETURN
149 END IF
150*
151* Quick return if possible
152*
153 IF( m.EQ.0 .OR. n.EQ.0 )
154 $ RETURN
155*
156* Determine the block size for this environment.
157*
158 nb = ilaenv( 1, 'DGETRF', ' ', m, n, -1, -1 )
159 IF( nb.LE.1 .OR. nb.GE.min( m, n ) ) THEN
160*
161* Use unblocked code.
162*
163 CALL dgetf2( m, n, a, lda, ipiv, info )
164 ELSE
165*
166* Use blocked code.
167*
168 DO 20 j = 1, min( m, n ), nb
169 jb = min( min( m, n )-j+1, nb )
170*
171* Update current block.
172*
173 CALL dgemm( 'No transpose', 'No transpose',
174 $ m-j+1, jb, j-1, -one,
175 $ a( j, 1 ), lda, a( 1, j ), lda, one,
176 $ a( j, j ), lda )
177
178*
179* Factor diagonal and subdiagonal blocks and test for exact
180* singularity.
181*
182 CALL dgetf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo )
183*
184* Adjust INFO and the pivot indices.
185*
186 IF( info.EQ.0 .AND. iinfo.GT.0 )
187 $ info = iinfo + j - 1
188 DO 10 i = j, min( m, j+jb-1 )
189 ipiv( i ) = j - 1 + ipiv( i )
190 10 CONTINUE
191*
192* Apply interchanges to column 1:J-1
193*
194 CALL dlaswp( j-1, a, lda, j, j+jb-1, ipiv, 1 )
195*
196 IF ( j+jb.LE.n ) THEN
197*
198* Apply interchanges to column J+JB:N
199*
200 CALL dlaswp( n-j-jb+1, a( 1, j+jb ), lda, j, j+jb-1,
201 $ ipiv, 1 )
202*
203 CALL dgemm( 'No transpose', 'No transpose',
204 $ jb, n-j-jb+1, j-1, -one,
205 $ a( j, 1 ), lda, a( 1, j+jb ), lda, one,
206 $ a( j, j+jb ), lda )
207*
208* Compute block row of U.
209*
210 CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit',
211 $ jb, n-j-jb+1, one, a( j, j ), lda,
212 $ a( j, j+jb ), lda )
213 END IF
214
215 20 CONTINUE
216
217 END IF
218 RETURN
219*
220* End of DGETRF
221*
dgetf2
subroutine dgetf2(m, n, a, lda, ipiv, info)
DGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row inter...
Definition dgetf2.f:108
dlaswp
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
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181

◆ sgeqrf()

subroutine sgeqrf ( integer m,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) tau,
real, dimension( * ) work,
integer lwork,
integer info )

SGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.

Purpose: 

!>
!> SGEQRF computes a QR factorization of a real M-by-N matrix A:
!> A = Q * R.
!>
!> This is the left-looking Level 3 BLAS version of the algorithm.
!>
!>
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 A.
!>          On exit, the elements on and above the diagonal of the array
!>          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
!>          upper triangular if m >= n); the elements below the diagonal,
!>          with the array TAU, represent the orthogonal matrix Q as a
!>          product of min(m,n) elementary reflectors (see Further
!>          Details).
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[out]TAU
!>          TAU is REAL array, dimension (min(M,N))
!>          The scalar factors of the elementary reflectors (see Further
!>          Details).
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!> 
!>          The dimension of the array WORK. The dimension can be divided into three parts.
!> 
!>          1) The part for the triangular factor T. If the very last T is not bigger
!>             than any of the rest, then this part is NB x ceiling(K/NB), otherwise,
!>             NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T
!> 
!>          2) The part for the very last T when T is bigger than any of the rest T.
!>             The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB,
!>             where K = min(M,N), NX is calculated by
!>                   NX = MAX( 0, ILAENV( 3, 'SGEQRF', ' ', M, N, -1, -1 ) )
!> 
!>          3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB)
!> 
!>          So LWORK = part1 + part2 + part3
!> 
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[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
December 2016

Further Details 

!>
!>  The matrix Q is represented as a product of elementary reflectors
!>
!>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
!>
!>  Each H(i) has the form
!>
!>     H(i) = I - tau * v * v'
!>
!>  where tau is a real scalar, and v is a real vector with
!>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
!>  and tau in TAU(i).
!>
!>

Definition at line 150 of file sgeqrf.f.

151*
152* -- LAPACK computational routine --
153* -- LAPACK is a software package provided by Univ. of Tennessee, --
154* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155*
156* .. Scalar Arguments ..
157 INTEGER INFO, LDA, LWORK, M, N
158* ..
159* .. Array Arguments ..
160 REAL A( LDA, * ), TAU( * ), WORK( * )
161* ..
162*
163* =====================================================================
164*
165* .. Local Scalars ..
166 LOGICAL LQUERY
167 INTEGER I, IB, IINFO, IWS, J, K, LWKOPT, NB,
168 $ NBMIN, NX, LBWORK, NT, LLWORK
169* ..
170* .. External Subroutines ..
171 EXTERNAL sgeqr2, slarfb, slarft, xerbla
172* ..
173* .. Intrinsic Functions ..
174 INTRINSIC max, min
175* ..
176* .. External Functions ..
177 INTEGER ILAENV
178 REAL SCEIL
179 EXTERNAL ilaenv, sceil
180* ..
181* .. Executable Statements ..
182
183 info = 0
184 nbmin = 2
185 nx = 0
186 iws = n
187 k = min( m, n )
188 nb = ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
189
190 IF( nb.GT.1 .AND. nb.LT.k ) THEN
191*
192* Determine when to cross over from blocked to unblocked code.
193*
194 nx = max( 0, ilaenv( 3, 'SGEQRF', ' ', m, n, -1, -1 ) )
195 END IF
196*
197* Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.:
198*
199* NB=3 2NB=6 K=10
200* | | |
201* 1--2--3--4--5--6--7--8--9--10
202* | \________/
203* K-NX=5 NT=4
204*
205* So here 4 x 4 is the last T stored in the workspace
206*
207 nt = k-sceil(real(k-nx)/real(nb))*nb
208
209*
210* optimal workspace = space for dlarfb + space for normal T's + space for the last T
211*
212 llwork = max(max((n-m)*k, (n-m)*nb), max(k*nb, nb*nb))
213 llwork = sceil(real(llwork)/real(nb))
214
215 IF ( nt.GT.nb ) THEN
216
217 lbwork = k-nt
218*
219* Optimal workspace for dlarfb = MAX(1,N)*NT
220*
221 lwkopt = (lbwork+llwork)*nb
222 work( 1 ) = (lwkopt+nt*nt)
223
224 ELSE
225
226 lbwork = sceil(real(k)/real(nb))*nb
227 lwkopt = (lbwork+llwork-nb)*nb
228 work( 1 ) = lwkopt
229
230 END IF
231
232*
233* Test the input arguments
234*
235 lquery = ( lwork.EQ.-1 )
236 IF( m.LT.0 ) THEN
237 info = -1
238 ELSE IF( n.LT.0 ) THEN
239 info = -2
240 ELSE IF( lda.LT.max( 1, m ) ) THEN
241 info = -4
242 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
243 info = -7
244 END IF
245 IF( info.NE.0 ) THEN
246 CALL xerbla( 'SGEQRF', -info )
247 RETURN
248 ELSE IF( lquery ) THEN
249 RETURN
250 END IF
251*
252* Quick return if possible
253*
254 IF( k.EQ.0 ) THEN
255 work( 1 ) = 1
256 RETURN
257 END IF
258*
259 IF( nb.GT.1 .AND. nb.LT.k ) THEN
260
261 IF( nx.LT.k ) THEN
262*
263* Determine if workspace is large enough for blocked code.
264*
265 IF ( nt.LE.nb ) THEN
266 iws = (lbwork+llwork-nb)*nb
267 ELSE
268 iws = (lbwork+llwork)*nb+nt*nt
269 END IF
270
271 IF( lwork.LT.iws ) THEN
272*
273* Not enough workspace to use optimal NB: reduce NB and
274* determine the minimum value of NB.
275*
276 IF ( nt.LE.nb ) THEN
277 nb = lwork / (llwork+(lbwork-nb))
278 ELSE
279 nb = (lwork-nt*nt)/(lbwork+llwork)
280 END IF
281
282 nbmin = max( 2, ilaenv( 2, 'SGEQRF', ' ', m, n, -1,
283 $ -1 ) )
284 END IF
285 END IF
286 END IF
287*
288 IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
289*
290* Use blocked code initially
291*
292 DO 10 i = 1, k - nx, nb
293 ib = min( k-i+1, nb )
294*
295* Update the current column using old T's
296*
297 DO 20 j = 1, i - nb, nb
298*
299* Apply H' to A(J:M,I:I+IB-1) from the left
300*
301 CALL slarfb( 'Left', 'Transpose', 'Forward',
302 $ 'Columnwise', m-j+1, ib, nb,
303 $ a( j, j ), lda, work(j), lbwork,
304 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
305 $ ib)
306
30720 CONTINUE
308*
309* Compute the QR factorization of the current block
310* A(I:M,I:I+IB-1)
311*
312 CALL sgeqr2( m-i+1, ib, a( i, i ), lda, tau( i ),
313 $ work(lbwork*nb+nt*nt+1), iinfo )
314
315 IF( i+ib.LE.n ) THEN
316*
317* Form the triangular factor of the block reflector
318* H = H(i) H(i+1) . . . H(i+ib-1)
319*
320 CALL slarft( 'Forward', 'Columnwise', m-i+1, ib,
321 $ a( i, i ), lda, tau( i ),
322 $ work(i), lbwork )
323*
324 END IF
325 10 CONTINUE
326 ELSE
327 i = 1
328 END IF
329*
330* Use unblocked code to factor the last or only block.
331*
332 IF( i.LE.k ) THEN
333
334 IF ( i .NE. 1 ) THEN
335
336 DO 30 j = 1, i - nb, nb
337*
338* Apply H' to A(J:M,I:K) from the left
339*
340 CALL slarfb( 'Left', 'Transpose', 'Forward',
341 $ 'Columnwise', m-j+1, k-i+1, nb,
342 $ a( j, j ), lda, work(j), lbwork,
343 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
344 $ k-i+1)
34530 CONTINUE
346
347 CALL sgeqr2( m-i+1, k-i+1, a( i, i ), lda, tau( i ),
348 $ work(lbwork*nb+nt*nt+1),iinfo )
349
350 ELSE
351*
352* Use unblocked code to factor the last or only block.
353*
354 CALL sgeqr2( m-i+1, n-i+1, a( i, i ), lda, tau( i ),
355 $ work,iinfo )
356
357 END IF
358 END IF
359
360
361*
362* Apply update to the column M+1:N when N > M
363*
364 IF ( m.LT.n .AND. i.NE.1) THEN
365*
366* Form the last triangular factor of the block reflector
367* H = H(i) H(i+1) . . . H(i+ib-1)
368*
369 IF ( nt .LE. nb ) THEN
370 CALL slarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
371 $ a( i, i ), lda, tau( i ), work(i), lbwork )
372 ELSE
373 CALL slarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
374 $ a( i, i ), lda, tau( i ),
375 $ work(lbwork*nb+1), nt )
376 END IF
377
378*
379* Apply H' to A(1:M,M+1:N) from the left
380*
381 DO 40 j = 1, k-nx, nb
382
383 ib = min( k-j+1, nb )
384
385 CALL slarfb( 'Left', 'Transpose', 'Forward',
386 $ 'Columnwise', m-j+1, n-m, ib,
387 $ a( j, j ), lda, work(j), lbwork,
388 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
389 $ n-m)
390
39140 CONTINUE
392
393 IF ( nt.LE.nb ) THEN
394 CALL slarfb( 'Left', 'Transpose', 'Forward',
395 $ 'Columnwise', m-j+1, n-m, k-j+1,
396 $ a( j, j ), lda, work(j), lbwork,
397 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
398 $ n-m)
399 ELSE
400 CALL slarfb( 'Left', 'Transpose', 'Forward',
401 $ 'Columnwise', m-j+1, n-m, k-j+1,
402 $ a( j, j ), lda,
403 $ work(lbwork*nb+1),
404 $ nt, a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
405 $ n-m)
406 END IF
407
408 END IF
409
410 work( 1 ) = iws
411 RETURN
412*
413* End of SGEQRF
414*
sgeqr2
subroutine sgeqr2(m, n, a, lda, tau, work, info)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition sgeqr2.f:130
slarft
subroutine slarft(direct, storev, n, k, v, ldv, tau, t, ldt)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition slarft.f:163
slarfb
subroutine slarfb(side, trans, direct, storev, m, n, k, v, ldv, t, ldt, c, ldc, work, ldwork)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition slarfb.f:197

◆ sgetrf()

subroutine sgetrf ( integer m,
integer n,
real, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer info )

SGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.

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

SGETRF VARIANT: left-looking Level 3 BLAS version of the 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 is the Crout Level 3 BLAS version of the algorithm.
!>
!>
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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

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 is the left-looking Level 3 BLAS version of the algorithm.
!>
!>
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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

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.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 101 of file sgetrf.f.

102*
103* -- LAPACK computational routine --
104* -- LAPACK is a software package provided by Univ. of Tennessee, --
105* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
106*
107* .. Scalar Arguments ..
108 INTEGER INFO, LDA, M, N
109* ..
110* .. Array Arguments ..
111 INTEGER IPIV( * )
112 REAL A( LDA, * )
113* ..
114*
115* =====================================================================
116*
117* .. Parameters ..
118 REAL ONE
119 parameter( one = 1.0e+0 )
120* ..
121* .. Local Scalars ..
122 INTEGER I, IINFO, J, JB, NB
123* ..
124* .. External Subroutines ..
125 EXTERNAL sgemm, sgetf2, slaswp, strsm, xerbla
126* ..
127* .. External Functions ..
128 INTEGER ILAENV
129 EXTERNAL ilaenv
130* ..
131* .. Intrinsic Functions ..
132 INTRINSIC max, min
133* ..
134* .. Executable Statements ..
135*
136* Test the input parameters.
137*
138 info = 0
139 IF( m.LT.0 ) THEN
140 info = -1
141 ELSE IF( n.LT.0 ) THEN
142 info = -2
143 ELSE IF( lda.LT.max( 1, m ) ) THEN
144 info = -4
145 END IF
146 IF( info.NE.0 ) THEN
147 CALL xerbla( 'SGETRF', -info )
148 RETURN
149 END IF
150*
151* Quick return if possible
152*
153 IF( m.EQ.0 .OR. n.EQ.0 )
154 $ RETURN
155*
156* Determine the block size for this environment.
157*
158 nb = ilaenv( 1, 'SGETRF', ' ', m, n, -1, -1 )
159 IF( nb.LE.1 .OR. nb.GE.min( m, n ) ) THEN
160*
161* Use unblocked code.
162*
163 CALL sgetf2( m, n, a, lda, ipiv, info )
164 ELSE
165*
166* Use blocked code.
167*
168 DO 20 j = 1, min( m, n ), nb
169 jb = min( min( m, n )-j+1, nb )
170*
171* Update current block.
172*
173 CALL sgemm( 'No transpose', 'No transpose',
174 $ m-j+1, jb, j-1, -one,
175 $ a( j, 1 ), lda, a( 1, j ), lda, one,
176 $ a( j, j ), lda )
177
178*
179* Factor diagonal and subdiagonal blocks and test for exact
180* singularity.
181*
182 CALL sgetf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo )
183*
184* Adjust INFO and the pivot indices.
185*
186 IF( info.EQ.0 .AND. iinfo.GT.0 )
187 $ info = iinfo + j - 1
188 DO 10 i = j, min( m, j+jb-1 )
189 ipiv( i ) = j - 1 + ipiv( i )
190 10 CONTINUE
191*
192* Apply interchanges to column 1:J-1
193*
194 CALL slaswp( j-1, a, lda, j, j+jb-1, ipiv, 1 )
195*
196 IF ( j+jb.LE.n ) THEN
197*
198* Apply interchanges to column J+JB:N
199*
200 CALL slaswp( n-j-jb+1, a( 1, j+jb ), lda, j, j+jb-1,
201 $ ipiv, 1 )
202*
203 CALL sgemm( 'No transpose', 'No transpose',
204 $ jb, n-j-jb+1, j-1, -one,
205 $ a( j, 1 ), lda, a( 1, j+jb ), lda, one,
206 $ a( j, j+jb ), lda )
207*
208* Compute block row of U.
209*
210 CALL strsm( 'Left', 'Lower', 'No transpose', 'Unit',
211 $ jb, n-j-jb+1, one, a( j, j ), lda,
212 $ a( j, j+jb ), lda )
213 END IF
214
215 20 CONTINUE
216
217 END IF
218 RETURN
219*
220* End of SGETRF
221*
sgetf2
subroutine sgetf2(m, n, a, lda, ipiv, info)
SGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row inter...
Definition sgetf2.f:108
slaswp
subroutine slaswp(n, a, lda, k1, k2, ipiv, incx)
SLASWP performs a series of row interchanges on a general rectangular matrix.
Definition slaswp.f:115
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187

◆ zgeqrf()

subroutine zgeqrf ( integer m,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( * ) tau,
complex*16, dimension( * ) work,
integer lwork,
integer info )

ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.

Purpose: 

!>
!> ZGEQRF computes a QR factorization of a real M-by-N matrix A:
!> A = Q * R.
!>
!> This is the left-looking Level 3 BLAS version of the algorithm.
!>
!>
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 A.
!>          On exit, the elements on and above the diagonal of the array
!>          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
!>          upper triangular if m >= n); the elements below the diagonal,
!>          with the array TAU, represent the orthogonal matrix Q as a
!>          product of min(m,n) elementary reflectors (see Further
!>          Details).
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[out]TAU
!>          TAU is COMPLEX*16 array, dimension (min(M,N))
!>          The scalar factors of the elementary reflectors (see Further
!>          Details).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!> 
!>          The dimension of the array WORK. The dimension can be divided into three parts.
!> 
!>          1) The part for the triangular factor T. If the very last T is not bigger
!>             than any of the rest, then this part is NB x ceiling(K/NB), otherwise,
!>             NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T
!> 
!>          2) The part for the very last T when T is bigger than any of the rest T.
!>             The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB,
!>             where K = min(M,N), NX is calculated by
!>                   NX = MAX( 0, ILAENV( 3, 'ZGEQRF', ' ', M, N, -1, -1 ) )
!> 
!>          3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB)
!> 
!>          So LWORK = part1 + part2 + part3
!> 
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[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
December 2016

Further Details 

!>
!>  The matrix Q is represented as a product of elementary reflectors
!>
!>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
!>
!>  Each H(i) has the form
!>
!>     H(i) = I - tau * v * v'
!>
!>  where tau is a real scalar, and v is a real vector with
!>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
!>  and tau in TAU(i).
!>
!>

Definition at line 150 of file zgeqrf.f.

151*
152* -- LAPACK computational routine --
153* -- LAPACK is a software package provided by Univ. of Tennessee, --
154* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155*
156* .. Scalar Arguments ..
157 INTEGER INFO, LDA, LWORK, M, N
158* ..
159* .. Array Arguments ..
160 COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
161* ..
162*
163* =====================================================================
164*
165* .. Local Scalars ..
166 LOGICAL LQUERY
167 INTEGER I, IB, IINFO, IWS, J, K, LWKOPT, NB,
168 $ NBMIN, NX, LBWORK, NT, LLWORK
169* ..
170* .. External Subroutines ..
171 EXTERNAL zgeqr2, zlarfb, zlarft, xerbla
172* ..
173* .. Intrinsic Functions ..
174 INTRINSIC max, min
175* ..
176* .. External Functions ..
177 INTEGER ILAENV
178 REAL SCEIL
179 EXTERNAL ilaenv, sceil
180* ..
181* .. Executable Statements ..
182
183 info = 0
184 nbmin = 2
185 nx = 0
186 iws = n
187 k = min( m, n )
188 nb = ilaenv( 1, 'ZGEQRF', ' ', m, n, -1, -1 )
189
190 IF( nb.GT.1 .AND. nb.LT.k ) THEN
191*
192* Determine when to cross over from blocked to unblocked code.
193*
194 nx = max( 0, ilaenv( 3, 'ZGEQRF', ' ', m, n, -1, -1 ) )
195 END IF
196*
197* Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.:
198*
199* NB=3 2NB=6 K=10
200* | | |
201* 1--2--3--4--5--6--7--8--9--10
202* | \________/
203* K-NX=5 NT=4
204*
205* So here 4 x 4 is the last T stored in the workspace
206*
207 nt = k-sceil(real(k-nx)/real(nb))*nb
208
209*
210* optimal workspace = space for dlarfb + space for normal T's + space for the last T
211*
212 llwork = max(max((n-m)*k, (n-m)*nb), max(k*nb, nb*nb))
213 llwork = sceil(real(llwork)/real(nb))
214
215 IF ( nt.GT.nb ) THEN
216
217 lbwork = k-nt
218*
219* Optimal workspace for dlarfb = MAX(1,N)*NT
220*
221 lwkopt = (lbwork+llwork)*nb
222 work( 1 ) = (lwkopt+nt*nt)
223
224 ELSE
225
226 lbwork = sceil(real(k)/real(nb))*nb
227 lwkopt = (lbwork+llwork-nb)*nb
228 work( 1 ) = lwkopt
229
230 END IF
231
232*
233* Test the input arguments
234*
235 lquery = ( lwork.EQ.-1 )
236 IF( m.LT.0 ) THEN
237 info = -1
238 ELSE IF( n.LT.0 ) THEN
239 info = -2
240 ELSE IF( lda.LT.max( 1, m ) ) THEN
241 info = -4
242 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
243 info = -7
244 END IF
245 IF( info.NE.0 ) THEN
246 CALL xerbla( 'ZGEQRF', -info )
247 RETURN
248 ELSE IF( lquery ) THEN
249 RETURN
250 END IF
251*
252* Quick return if possible
253*
254 IF( k.EQ.0 ) THEN
255 work( 1 ) = 1
256 RETURN
257 END IF
258*
259 IF( nb.GT.1 .AND. nb.LT.k ) THEN
260
261 IF( nx.LT.k ) THEN
262*
263* Determine if workspace is large enough for blocked code.
264*
265 IF ( nt.LE.nb ) THEN
266 iws = (lbwork+llwork-nb)*nb
267 ELSE
268 iws = (lbwork+llwork)*nb+nt*nt
269 END IF
270
271 IF( lwork.LT.iws ) THEN
272*
273* Not enough workspace to use optimal NB: reduce NB and
274* determine the minimum value of NB.
275*
276 IF ( nt.LE.nb ) THEN
277 nb = lwork / (llwork+(lbwork-nb))
278 ELSE
279 nb = (lwork-nt*nt)/(lbwork+llwork)
280 END IF
281
282 nbmin = max( 2, ilaenv( 2, 'ZGEQRF', ' ', m, n, -1,
283 $ -1 ) )
284 END IF
285 END IF
286 END IF
287*
288 IF( nb.GE.nbmin .AND. nb.LT.k .AND. nx.LT.k ) THEN
289*
290* Use blocked code initially
291*
292 DO 10 i = 1, k - nx, nb
293 ib = min( k-i+1, nb )
294*
295* Update the current column using old T's
296*
297 DO 20 j = 1, i - nb, nb
298*
299* Apply H' to A(J:M,I:I+IB-1) from the left
300*
301 CALL zlarfb( 'Left', 'Transpose', 'Forward',
302 $ 'Columnwise', m-j+1, ib, nb,
303 $ a( j, j ), lda, work(j), lbwork,
304 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
305 $ ib)
306
30720 CONTINUE
308*
309* Compute the QR factorization of the current block
310* A(I:M,I:I+IB-1)
311*
312 CALL zgeqr2( m-i+1, ib, a( i, i ), lda, tau( i ),
313 $ work(lbwork*nb+nt*nt+1), iinfo )
314
315 IF( i+ib.LE.n ) THEN
316*
317* Form the triangular factor of the block reflector
318* H = H(i) H(i+1) . . . H(i+ib-1)
319*
320 CALL zlarft( 'Forward', 'Columnwise', m-i+1, ib,
321 $ a( i, i ), lda, tau( i ),
322 $ work(i), lbwork )
323*
324 END IF
325 10 CONTINUE
326 ELSE
327 i = 1
328 END IF
329*
330* Use unblocked code to factor the last or only block.
331*
332 IF( i.LE.k ) THEN
333
334 IF ( i .NE. 1 ) THEN
335
336 DO 30 j = 1, i - nb, nb
337*
338* Apply H' to A(J:M,I:K) from the left
339*
340 CALL zlarfb( 'Left', 'Transpose', 'Forward',
341 $ 'Columnwise', m-j+1, k-i+1, nb,
342 $ a( j, j ), lda, work(j), lbwork,
343 $ a( j, i ), lda, work(lbwork*nb+nt*nt+1),
344 $ k-i+1)
34530 CONTINUE
346
347 CALL zgeqr2( m-i+1, k-i+1, a( i, i ), lda, tau( i ),
348 $ work(lbwork*nb+nt*nt+1),iinfo )
349
350 ELSE
351*
352* Use unblocked code to factor the last or only block.
353*
354 CALL zgeqr2( m-i+1, n-i+1, a( i, i ), lda, tau( i ),
355 $ work,iinfo )
356
357 END IF
358 END IF
359
360
361*
362* Apply update to the column M+1:N when N > M
363*
364 IF ( m.LT.n .AND. i.NE.1) THEN
365*
366* Form the last triangular factor of the block reflector
367* H = H(i) H(i+1) . . . H(i+ib-1)
368*
369 IF ( nt .LE. nb ) THEN
370 CALL zlarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
371 $ a( i, i ), lda, tau( i ), work(i), lbwork )
372 ELSE
373 CALL zlarft( 'Forward', 'Columnwise', m-i+1, k-i+1,
374 $ a( i, i ), lda, tau( i ),
375 $ work(lbwork*nb+1), nt )
376 END IF
377
378*
379* Apply H' to A(1:M,M+1:N) from the left
380*
381 DO 40 j = 1, k-nx, nb
382
383 ib = min( k-j+1, nb )
384
385 CALL zlarfb( 'Left', 'Transpose', 'Forward',
386 $ 'Columnwise', m-j+1, n-m, ib,
387 $ a( j, j ), lda, work(j), lbwork,
388 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
389 $ n-m)
390
39140 CONTINUE
392
393 IF ( nt.LE.nb ) THEN
394 CALL zlarfb( 'Left', 'Transpose', 'Forward',
395 $ 'Columnwise', m-j+1, n-m, k-j+1,
396 $ a( j, j ), lda, work(j), lbwork,
397 $ a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
398 $ n-m)
399 ELSE
400 CALL zlarfb( 'Left', 'Transpose', 'Forward',
401 $ 'Columnwise', m-j+1, n-m, k-j+1,
402 $ a( j, j ), lda,
403 $ work(lbwork*nb+1),
404 $ nt, a( j, m+1 ), lda, work(lbwork*nb+nt*nt+1),
405 $ n-m)
406 END IF
407
408 END IF
409
410 work( 1 ) = iws
411 RETURN
412*
413* End of ZGEQRF
414*
zgeqr2
subroutine zgeqr2(m, n, a, lda, tau, work, info)
ZGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition zgeqr2.f:130
zlarfb
subroutine zlarfb(side, trans, direct, storev, m, n, k, v, ldv, t, ldt, c, ldc, work, ldwork)
ZLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
Definition zlarfb.f:197
zlarft
subroutine zlarft(direct, storev, n, k, v, ldv, tau, t, ldt)
ZLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition zlarft.f:163

◆ zgetrf()

subroutine zgetrf ( integer m,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer info )

ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.

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

ZGETRF VARIANT: left-looking Level 3 BLAS version of the 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 is the Crout Level 3 BLAS version of the algorithm.
!>
!>
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
December 2016

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 is the left-looking Level 3 BLAS version of the algorithm.
!>
!>
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
December 2016

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
December 2016

Definition at line 101 of file zgetrf.f.

102*
103* -- LAPACK computational routine --
104* -- LAPACK is a software package provided by Univ. of Tennessee, --
105* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
106*
107* .. Scalar Arguments ..
108 INTEGER INFO, LDA, M, N
109* ..
110* .. Array Arguments ..
111 INTEGER IPIV( * )
112 COMPLEX*16 A( LDA, * )
113* ..
114*
115* =====================================================================
116*
117* .. Parameters ..
118 COMPLEX*16 ONE
119 parameter( one = ( 1.0d+0, 0.0d+0 ) )
120* ..
121* .. Local Scalars ..
122 INTEGER I, IINFO, J, JB, NB
123* ..
124* .. External Subroutines ..
125 EXTERNAL zgemm, zgetf2, zlaswp, ztrsm, xerbla
126* ..
127* .. External Functions ..
128 INTEGER ILAENV
129 EXTERNAL ilaenv
130* ..
131* .. Intrinsic Functions ..
132 INTRINSIC max, min
133* ..
134* .. Executable Statements ..
135*
136* Test the input parameters.
137*
138 info = 0
139 IF( m.LT.0 ) THEN
140 info = -1
141 ELSE IF( n.LT.0 ) THEN
142 info = -2
143 ELSE IF( lda.LT.max( 1, m ) ) THEN
144 info = -4
145 END IF
146 IF( info.NE.0 ) THEN
147 CALL xerbla( 'ZGETRF', -info )
148 RETURN
149 END IF
150*
151* Quick return if possible
152*
153 IF( m.EQ.0 .OR. n.EQ.0 )
154 $ RETURN
155*
156* Determine the block size for this environment.
157*
158 nb = ilaenv( 1, 'ZGETRF', ' ', m, n, -1, -1 )
159 IF( nb.LE.1 .OR. nb.GE.min( m, n ) ) THEN
160*
161* Use unblocked code.
162*
163 CALL zgetf2( m, n, a, lda, ipiv, info )
164 ELSE
165*
166* Use blocked code.
167*
168 DO 20 j = 1, min( m, n ), nb
169 jb = min( min( m, n )-j+1, nb )
170*
171* Update current block.
172*
173 CALL zgemm( 'No transpose', 'No transpose',
174 $ m-j+1, jb, j-1, -one,
175 $ a( j, 1 ), lda, a( 1, j ), lda, one,
176 $ a( j, j ), lda )
177
178*
179* Factor diagonal and subdiagonal blocks and test for exact
180* singularity.
181*
182 CALL zgetf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo )
183*
184* Adjust INFO and the pivot indices.
185*
186 IF( info.EQ.0 .AND. iinfo.GT.0 )
187 $ info = iinfo + j - 1
188 DO 10 i = j, min( m, j+jb-1 )
189 ipiv( i ) = j - 1 + ipiv( i )
190 10 CONTINUE
191*
192* Apply interchanges to column 1:J-1
193*
194 CALL zlaswp( j-1, a, lda, j, j+jb-1, ipiv, 1 )
195*
196 IF ( j+jb.LE.n ) THEN
197*
198* Apply interchanges to column J+JB:N
199*
200 CALL zlaswp( n-j-jb+1, a( 1, j+jb ), lda, j, j+jb-1,
201 $ ipiv, 1 )
202*
203 CALL zgemm( 'No transpose', 'No transpose',
204 $ jb, n-j-jb+1, j-1, -one,
205 $ a( j, 1 ), lda, a( 1, j+jb ), lda, one,
206 $ a( j, j+jb ), lda )
207*
208* Compute block row of U.
209*
210 CALL ztrsm( 'Left', 'Lower', 'No transpose', 'Unit',
211 $ jb, n-j-jb+1, one, a( j, j ), lda,
212 $ a( j, j+jb ), lda )
213 END IF
214
215 20 CONTINUE
216
217 END IF
218 RETURN
219*
220* End of ZGETRF
221*
zgetf2
subroutine zgetf2(m, n, a, lda, ipiv, info)
ZGETF2 computes the LU factorization of a general m-by-n matrix using partial pivoting with row inter...
Definition zgetf2.f:108
zlaswp
subroutine zlaswp(n, a, lda, k1, k2, ipiv, incx)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
Definition zlaswp.f:115
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:187
ztrsm
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180