145 SUBROUTINE dgbtrf( M, N, KL, KU, AB, LDAB, IPIV, INFO )
153 INTEGER info, kl, ku, ldab, m, n
157 DOUBLE PRECISION ab( ldab, * )
163 DOUBLE PRECISION one, zero
164 parameter( one = 1.0d+0, zero = 0.0d+0 )
165 INTEGER nbmax, ldwork
166 parameter( nbmax = 64, ldwork = nbmax+1 )
169 INTEGER i, i2, i3, ii, ip,
j, j2, j3, jb, jj, jm, jp,
170 $ ju, k2, km, kv, nb, nw
171 DOUBLE PRECISION temp
174 DOUBLE PRECISION work13( ldwork, nbmax ),
175 $ work31( ldwork, nbmax )
200 ELSE IF( n.LT.0 )
THEN
202 ELSE IF( kl.LT.0 )
THEN
204 ELSE IF( ku.LT.0 )
THEN
206 ELSE IF( ldab.LT.kl+kv+1 )
THEN
210 CALL
xerbla(
'DGBTRF', -info )
216 IF( m.EQ.0 .OR. n.EQ.0 )
221 nb =
ilaenv( 1,
'DGBTRF',
' ', m, n, kl, ku )
226 nb = min( nb, nbmax )
228 IF( nb.LE.1 .OR. nb.GT.kl )
THEN
232 CALL
dgbtf2( m, n, kl, ku, ab, ldab, ipiv, info )
241 work13( i,
j ) = zero
249 work31( i,
j ) = zero
257 DO 60
j = ku + 2, min( kv, n )
258 DO 50 i = kv -
j + 2, kl
268 DO 180
j = 1, min( m, n ), nb
269 jb = min( nb, min( m, n )-
j+1 )
283 i2 = min( kl-jb, m-
j-jb+1 )
284 i3 = min( jb, m-
j-kl+1 )
290 DO 80 jj =
j,
j + jb - 1
294 IF( jj+kv.LE.n )
THEN
296 ab( i, jj+kv ) = zero
304 jp =
idamax( km+1, ab( kv+1, jj ), 1 )
305 ipiv( jj ) = jp + jj -
j
306 IF( ab( kv+jp, jj ).NE.zero )
THEN
307 ju = max( ju, min( jj+ku+jp-1, n ) )
312 IF( jp+jj-1.LT.
j+kl )
THEN
314 CALL
dswap( jb, ab( kv+1+jj-
j,
j ), ldab-1,
315 $ ab( kv+jp+jj-
j,
j ), ldab-1 )
321 CALL
dswap( jj-
j, ab( kv+1+jj-
j,
j ), ldab-1,
322 $ work31( jp+jj-
j-kl, 1 ), ldwork )
323 CALL
dswap(
j+jb-jj, ab( kv+1, jj ), ldab-1,
324 $ ab( kv+jp, jj ), ldab-1 )
330 CALL
dscal( km, one / ab( kv+1, jj ), ab( kv+2, jj ),
337 jm = min( ju,
j+jb-1 )
339 $ CALL
dger( km, jm-jj, -one, ab( kv+2, jj ), 1,
340 $ ab( kv, jj+1 ), ldab-1,
341 $ ab( kv+1, jj+1 ), ldab-1 )
353 nw = min( jj-
j+1, i3 )
355 $ CALL
dcopy( nw, ab( kv+kl+1-jj+
j, jj ), 1,
356 $ work31( 1, jj-
j+1 ), 1 )
362 j2 = min( ju-
j+1, kv ) - jb
363 j3 = max( 0, ju-
j-kv+1 )
368 CALL
dlaswp( j2, ab( kv+1-jb,
j+jb ), ldab-1, 1, jb,
373 DO 90 i =
j,
j + jb - 1
374 ipiv( i ) = ipiv( i ) +
j - 1
383 DO 100 ii =
j + i - 1,
j + jb - 1
386 temp = ab( kv+1+ii-jj, jj )
387 ab( kv+1+ii-jj, jj ) = ab( kv+1+ip-jj, jj )
388 ab( kv+1+ip-jj, jj ) = temp
399 CALL
dtrsm(
'Left',
'Lower',
'No transpose',
'Unit',
400 $ jb, j2, one, ab( kv+1,
j ), ldab-1,
401 $ ab( kv+1-jb,
j+jb ), ldab-1 )
407 CALL
dgemm(
'No transpose',
'No transpose', i2, j2,
408 $ jb, -one, ab( kv+1+jb,
j ), ldab-1,
409 $ ab( kv+1-jb,
j+jb ), ldab-1, one,
410 $ ab( kv+1,
j+jb ), ldab-1 )
417 CALL
dgemm(
'No transpose',
'No transpose', i3, j2,
418 $ jb, -one, work31, ldwork,
419 $ ab( kv+1-jb,
j+jb ), ldab-1, one,
420 $ ab( kv+kl+1-jb,
j+jb ), ldab-1 )
431 work13( ii, jj ) = ab( ii-jj+1, jj+
j+kv-1 )
437 CALL
dtrsm(
'Left',
'Lower',
'No transpose',
'Unit',
438 $ jb, j3, one, ab( kv+1,
j ), ldab-1,
445 CALL
dgemm(
'No transpose',
'No transpose', i2, j3,
446 $ jb, -one, ab( kv+1+jb,
j ), ldab-1,
447 $ work13, ldwork, one, ab( 1+jb,
j+kv ),
455 CALL
dgemm(
'No transpose',
'No transpose', i3, j3,
456 $ jb, -one, work31, ldwork, work13,
457 $ ldwork, one, ab( 1+kl,
j+kv ), ldab-1 )
464 ab( ii-jj+1, jj+
j+kv-1 ) = work13( ii, jj )
472 DO 160 i =
j,
j + jb - 1
473 ipiv( i ) = ipiv( i ) +
j - 1
481 DO 170 jj =
j + jb - 1,
j, -1
482 jp = ipiv( jj ) - jj + 1
487 IF( jp+jj-1.LT.
j+kl )
THEN
491 CALL
dswap( jj-
j, ab( kv+1+jj-
j,
j ), ldab-1,
492 $ ab( kv+jp+jj-
j,
j ), ldab-1 )
497 CALL
dswap( jj-
j, ab( kv+1+jj-
j,
j ), ldab-1,
498 $ work31( jp+jj-
j-kl, 1 ), ldwork )
504 nw = min( i3, jj-
j+1 )
506 $ CALL
dcopy( nw, work31( 1, jj-
j+1 ), 1,
507 $ ab( kv+kl+1-jj+
j, jj ), 1 )
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
DLASWP performs a series of row interchanges on a general rectangular matrix.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
DGBTRF
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
subroutine dgbtf2(M, N, KL, KU, AB, LDAB, IPIV, INFO)
DGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algo...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
integer function idamax(N, DX, INCX)
IDAMAX
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)