145 SUBROUTINE cgbtrf( M, N, KL, KU, AB, LDAB, IPIV, INFO )
153 INTEGER INFO, KL, KU, LDAB, M, N
157 COMPLEX AB( ldab, * )
164 parameter ( one = ( 1.0e+0, 0.0e+0 ),
165 $ zero = ( 0.0e+0, 0.0e+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 WORK13( ldwork, nbmax ),
176 $ work31( ldwork, nbmax )
179 INTEGER ICAMAX, ILAENV
180 EXTERNAL icamax, ilaenv
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(
'CGBTRF', -info )
217 IF( m.EQ.0 .OR. n.EQ.0 )
222 nb = ilaenv( 1,
'CGBTRF',
' ', m, n, kl, ku )
227 nb = min( nb, nbmax )
229 IF( nb.LE.1 .OR. nb.GT.kl )
THEN
233 CALL cgbtf2( 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 = icamax( 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 cswap( jb, ab( kv+1+jj-j, j ), ldab-1,
316 $ ab( kv+jp+jj-j, j ), ldab-1 )
322 CALL cswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
323 $ work31( jp+jj-j-kl, 1 ), ldwork )
324 CALL cswap( j+jb-jj, ab( kv+1, jj ), ldab-1,
325 $ ab( kv+jp, jj ), ldab-1 )
331 CALL cscal( km, one / ab( kv+1, jj ), ab( kv+2, jj ),
338 jm = min( ju, j+jb-1 )
340 $
CALL cgeru( 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 ccopy( 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 claswp( 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 ctrsm(
'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 cgemm(
'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 cgemm(
'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 ctrsm(
'Left',
'Lower',
'No transpose',
'Unit',
439 $ jb, j3, one, ab( kv+1, j ), ldab-1,
446 CALL cgemm(
'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 cgemm(
'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 cswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
493 $ ab( kv+jp+jj-j, j ), ldab-1 )
498 CALL cswap( 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 ccopy( nw, work31( 1, jj-j+1 ), 1,
508 $ ab( kv+kl+1-jj+j, jj ), 1 )
subroutine cgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
CGBTRF
subroutine claswp(N, A, LDA, K1, K2, IPIV, INCX)
CLASWP performs a series of row interchanges on a general rectangular matrix.
subroutine cscal(N, CA, CX, INCX)
CSCAL
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine cgbtf2(M, N, KL, KU, AB, LDAB, IPIV, INFO)
CGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algo...
subroutine cgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
CGERU
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM