145 SUBROUTINE zgbtrf( M, N, KL, KU, AB, LDAB, IPIV, INFO )
153 INTEGER INFO, KL, KU, LDAB, M, N
157 COMPLEX*16 AB( ldab, * )
164 parameter ( one = ( 1.0d+0, 0.0d+0 ),
165 $ zero = ( 0.0d+0, 0.0d+0 ) )
166 INTEGER NBMAX, LDWORK
167 parameter ( nbmax = 64, ldwork = nbmax+1 )
170 INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP,
171 $ ju, k2, km, kv, nb, nw
175 COMPLEX*16 WORK13( ldwork, nbmax ),
176 $ work31( ldwork, nbmax )
179 INTEGER ILAENV, IZAMAX
180 EXTERNAL ilaenv, izamax
201 ELSE IF( n.LT.0 )
THEN
203 ELSE IF( kl.LT.0 )
THEN
205 ELSE IF( ku.LT.0 )
THEN
207 ELSE IF( ldab.LT.kl+kv+1 )
THEN
211 CALL xerbla(
'ZGBTRF', -info )
217 IF( m.EQ.0 .OR. n.EQ.0 )
222 nb = ilaenv( 1,
'ZGBTRF',
' ', m, n, kl, ku )
227 nb = min( nb, nbmax )
229 IF( nb.LE.1 .OR. nb.GT.kl )
THEN
233 CALL zgbtf2( m, n, kl, ku, ab, ldab, ipiv, info )
242 work13( i, j ) = zero
250 work31( i, j ) = zero
258 DO 60 j = ku + 2, min( kv, n )
259 DO 50 i = kv - j + 2, kl
269 DO 180 j = 1, min( m, n ), nb
270 jb = min( nb, min( m, n )-j+1 )
284 i2 = min( kl-jb, m-j-jb+1 )
285 i3 = min( jb, m-j-kl+1 )
291 DO 80 jj = j, j + jb - 1
295 IF( jj+kv.LE.n )
THEN
297 ab( i, jj+kv ) = zero
305 jp = izamax( km+1, ab( kv+1, jj ), 1 )
306 ipiv( jj ) = jp + jj - j
307 IF( ab( kv+jp, jj ).NE.zero )
THEN
308 ju = max( ju, min( jj+ku+jp-1, n ) )
313 IF( jp+jj-1.LT.j+kl )
THEN
315 CALL zswap( jb, ab( kv+1+jj-j, j ), ldab-1,
316 $ ab( kv+jp+jj-j, j ), ldab-1 )
322 CALL zswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
323 $ work31( jp+jj-j-kl, 1 ), ldwork )
324 CALL zswap( j+jb-jj, ab( kv+1, jj ), ldab-1,
325 $ ab( kv+jp, jj ), ldab-1 )
331 CALL zscal( km, one / ab( kv+1, jj ), ab( kv+2, jj ),
338 jm = min( ju, j+jb-1 )
340 $
CALL zgeru( km, jm-jj, -one, ab( kv+2, jj ), 1,
341 $ ab( kv, jj+1 ), ldab-1,
342 $ ab( kv+1, jj+1 ), ldab-1 )
354 nw = min( jj-j+1, i3 )
356 $
CALL zcopy( nw, ab( kv+kl+1-jj+j, jj ), 1,
357 $ work31( 1, jj-j+1 ), 1 )
363 j2 = min( ju-j+1, kv ) - jb
364 j3 = max( 0, ju-j-kv+1 )
369 CALL zlaswp( j2, ab( kv+1-jb, j+jb ), ldab-1, 1, jb,
374 DO 90 i = j, j + jb - 1
375 ipiv( i ) = ipiv( i ) + j - 1
384 DO 100 ii = j + i - 1, j + jb - 1
387 temp = ab( kv+1+ii-jj, jj )
388 ab( kv+1+ii-jj, jj ) = ab( kv+1+ip-jj, jj )
389 ab( kv+1+ip-jj, jj ) = temp
400 CALL ztrsm(
'Left',
'Lower',
'No transpose',
'Unit',
401 $ jb, j2, one, ab( kv+1, j ), ldab-1,
402 $ ab( kv+1-jb, j+jb ), ldab-1 )
408 CALL zgemm(
'No transpose',
'No transpose', i2, j2,
409 $ jb, -one, ab( kv+1+jb, j ), ldab-1,
410 $ ab( kv+1-jb, j+jb ), ldab-1, one,
411 $ ab( kv+1, j+jb ), ldab-1 )
418 CALL zgemm(
'No transpose',
'No transpose', i3, j2,
419 $ jb, -one, work31, ldwork,
420 $ ab( kv+1-jb, j+jb ), ldab-1, one,
421 $ ab( kv+kl+1-jb, j+jb ), ldab-1 )
432 work13( ii, jj ) = ab( ii-jj+1, jj+j+kv-1 )
438 CALL ztrsm(
'Left',
'Lower',
'No transpose',
'Unit',
439 $ jb, j3, one, ab( kv+1, j ), ldab-1,
446 CALL zgemm(
'No transpose',
'No transpose', i2, j3,
447 $ jb, -one, ab( kv+1+jb, j ), ldab-1,
448 $ work13, ldwork, one, ab( 1+jb, j+kv ),
456 CALL zgemm(
'No transpose',
'No transpose', i3, j3,
457 $ jb, -one, work31, ldwork, work13,
458 $ ldwork, one, ab( 1+kl, j+kv ), ldab-1 )
465 ab( ii-jj+1, jj+j+kv-1 ) = work13( ii, jj )
473 DO 160 i = j, j + jb - 1
474 ipiv( i ) = ipiv( i ) + j - 1
482 DO 170 jj = j + jb - 1, j, -1
483 jp = ipiv( jj ) - jj + 1
488 IF( jp+jj-1.LT.j+kl )
THEN
492 CALL zswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
493 $ ab( kv+jp+jj-j, j ), ldab-1 )
498 CALL zswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
499 $ work31( jp+jj-j-kl, 1 ), ldwork )
505 nw = min( i3, jj-j+1 )
507 $
CALL zcopy( nw, work31( 1, jj-j+1 ), 1,
508 $ ab( kv+kl+1-jj+j, jj ), 1 )
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zlaswp(N, A, LDA, K1, K2, IPIV, INCX)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
ZGBTRF
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zgbtf2(M, N, KL, KU, AB, LDAB, IPIV, INFO)
ZGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algo...
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
subroutine zgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERU
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL