LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgbtrf()

subroutine dgbtrf ( integer  m,
integer  n,
integer  kl,
integer  ku,
double precision, dimension( ldab, * )  ab,
integer  ldab,
integer, dimension( * )  ipiv,
integer  info 
)

DGBTRF

Download DGBTRF + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DGBTRF computes an LU factorization of a real m-by-n band matrix A
 using partial pivoting with row interchanges.

 This is the blocked version of the algorithm, calling Level 3 BLAS.
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]KL
          KL is INTEGER
          The number of subdiagonals within the band of A.  KL >= 0.
[in]KU
          KU is INTEGER
          The number of superdiagonals within the band of A.  KU >= 0.
[in,out]AB
          AB is DOUBLE PRECISION array, dimension (LDAB,N)
          On entry, the matrix A in band storage, in rows KL+1 to
          2*KL+KU+1; rows 1 to KL of the array need not be set.
          The j-th column of A is stored in the j-th column of the
          array AB as follows:
          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)

          On exit, details of the factorization: U is stored as an
          upper triangular band matrix with KL+KU superdiagonals in
          rows 1 to KL+KU+1, and the multipliers used during the
          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
          See below for further details.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
[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.
Further Details:
  The band storage scheme is illustrated by the following example, when
  M = N = 6, KL = 2, KU = 1:

  On entry:                       On exit:

      *    *    *    +    +    +       *    *    *   u14  u25  u36
      *    *    +    +    +    +       *    *   u13  u24  u35  u46
      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *

  Array elements marked * are not used by the routine; elements marked
  + need not be set on entry, but are required by the routine to store
  elements of U because of fill-in resulting from the row interchanges.

Definition at line 143 of file dgbtrf.f.

144*
145* -- LAPACK computational routine --
146* -- LAPACK is a software package provided by Univ. of Tennessee, --
147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148*
149* .. Scalar Arguments ..
150 INTEGER INFO, KL, KU, LDAB, M, N
151* ..
152* .. Array Arguments ..
153 INTEGER IPIV( * )
154 DOUBLE PRECISION AB( LDAB, * )
155* ..
156*
157* =====================================================================
158*
159* .. Parameters ..
160 DOUBLE PRECISION ONE, ZERO
161 parameter( one = 1.0d+0, zero = 0.0d+0 )
162 INTEGER NBMAX, LDWORK
163 parameter( nbmax = 64, ldwork = nbmax+1 )
164* ..
165* .. Local Scalars ..
166 INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP,
167 $ JU, K2, KM, KV, NB, NW
168 DOUBLE PRECISION TEMP
169* ..
170* .. Local Arrays ..
171 DOUBLE PRECISION WORK13( LDWORK, NBMAX ),
172 $ WORK31( LDWORK, NBMAX )
173* ..
174* .. External Functions ..
175 INTEGER IDAMAX, ILAENV
176 EXTERNAL idamax, ilaenv
177* ..
178* .. External Subroutines ..
179 EXTERNAL dcopy, dgbtf2, dgemm, dger, dlaswp, dscal,
180 $ dswap, dtrsm, xerbla
181* ..
182* .. Intrinsic Functions ..
183 INTRINSIC max, min
184* ..
185* .. Executable Statements ..
186*
187* KV is the number of superdiagonals in the factor U, allowing for
188* fill-in
189*
190 kv = ku + kl
191*
192* Test the input parameters.
193*
194 info = 0
195 IF( m.LT.0 ) THEN
196 info = -1
197 ELSE IF( n.LT.0 ) THEN
198 info = -2
199 ELSE IF( kl.LT.0 ) THEN
200 info = -3
201 ELSE IF( ku.LT.0 ) THEN
202 info = -4
203 ELSE IF( ldab.LT.kl+kv+1 ) THEN
204 info = -6
205 END IF
206 IF( info.NE.0 ) THEN
207 CALL xerbla( 'DGBTRF', -info )
208 RETURN
209 END IF
210*
211* Quick return if possible
212*
213 IF( m.EQ.0 .OR. n.EQ.0 )
214 $ RETURN
215*
216* Determine the block size for this environment
217*
218 nb = ilaenv( 1, 'DGBTRF', ' ', m, n, kl, ku )
219*
220* The block size must not exceed the limit set by the size of the
221* local arrays WORK13 and WORK31.
222*
223 nb = min( nb, nbmax )
224*
225 IF( nb.LE.1 .OR. nb.GT.kl ) THEN
226*
227* Use unblocked code
228*
229 CALL dgbtf2( m, n, kl, ku, ab, ldab, ipiv, info )
230 ELSE
231*
232* Use blocked code
233*
234* Zero the superdiagonal elements of the work array WORK13
235*
236 DO 20 j = 1, nb
237 DO 10 i = 1, j - 1
238 work13( i, j ) = zero
239 10 CONTINUE
240 20 CONTINUE
241*
242* Zero the subdiagonal elements of the work array WORK31
243*
244 DO 40 j = 1, nb
245 DO 30 i = j + 1, nb
246 work31( i, j ) = zero
247 30 CONTINUE
248 40 CONTINUE
249*
250* Gaussian elimination with partial pivoting
251*
252* Set fill-in elements in columns KU+2 to KV to zero
253*
254 DO 60 j = ku + 2, min( kv, n )
255 DO 50 i = kv - j + 2, kl
256 ab( i, j ) = zero
257 50 CONTINUE
258 60 CONTINUE
259*
260* JU is the index of the last column affected by the current
261* stage of the factorization
262*
263 ju = 1
264*
265 DO 180 j = 1, min( m, n ), nb
266 jb = min( nb, min( m, n )-j+1 )
267*
268* The active part of the matrix is partitioned
269*
270* A11 A12 A13
271* A21 A22 A23
272* A31 A32 A33
273*
274* Here A11, A21 and A31 denote the current block of JB columns
275* which is about to be factorized. The number of rows in the
276* partitioning are JB, I2, I3 respectively, and the numbers
277* of columns are JB, J2, J3. The superdiagonal elements of A13
278* and the subdiagonal elements of A31 lie outside the band.
279*
280 i2 = min( kl-jb, m-j-jb+1 )
281 i3 = min( jb, m-j-kl+1 )
282*
283* J2 and J3 are computed after JU has been updated.
284*
285* Factorize the current block of JB columns
286*
287 DO 80 jj = j, j + jb - 1
288*
289* Set fill-in elements in column JJ+KV to zero
290*
291 IF( jj+kv.LE.n ) THEN
292 DO 70 i = 1, kl
293 ab( i, jj+kv ) = zero
294 70 CONTINUE
295 END IF
296*
297* Find pivot and test for singularity. KM is the number of
298* subdiagonal elements in the current column.
299*
300 km = min( kl, m-jj )
301 jp = idamax( km+1, ab( kv+1, jj ), 1 )
302 ipiv( jj ) = jp + jj - j
303 IF( ab( kv+jp, jj ).NE.zero ) THEN
304 ju = max( ju, min( jj+ku+jp-1, n ) )
305 IF( jp.NE.1 ) THEN
306*
307* Apply interchange to columns J to J+JB-1
308*
309 IF( jp+jj-1.LT.j+kl ) THEN
310*
311 CALL dswap( jb, ab( kv+1+jj-j, j ), ldab-1,
312 $ ab( kv+jp+jj-j, j ), ldab-1 )
313 ELSE
314*
315* The interchange affects columns J to JJ-1 of A31
316* which are stored in the work array WORK31
317*
318 CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
319 $ work31( jp+jj-j-kl, 1 ), ldwork )
320 CALL dswap( j+jb-jj, ab( kv+1, jj ), ldab-1,
321 $ ab( kv+jp, jj ), ldab-1 )
322 END IF
323 END IF
324*
325* Compute multipliers
326*
327 CALL dscal( km, one / ab( kv+1, jj ), ab( kv+2, jj ),
328 $ 1 )
329*
330* Update trailing submatrix within the band and within
331* the current block. JM is the index of the last column
332* which needs to be updated.
333*
334 jm = min( ju, j+jb-1 )
335 IF( jm.GT.jj )
336 $ CALL dger( km, jm-jj, -one, ab( kv+2, jj ), 1,
337 $ ab( kv, jj+1 ), ldab-1,
338 $ ab( kv+1, jj+1 ), ldab-1 )
339 ELSE
340*
341* If pivot is zero, set INFO to the index of the pivot
342* unless a zero pivot has already been found.
343*
344 IF( info.EQ.0 )
345 $ info = jj
346 END IF
347*
348* Copy current column of A31 into the work array WORK31
349*
350 nw = min( jj-j+1, i3 )
351 IF( nw.GT.0 )
352 $ CALL dcopy( nw, ab( kv+kl+1-jj+j, jj ), 1,
353 $ work31( 1, jj-j+1 ), 1 )
354 80 CONTINUE
355 IF( j+jb.LE.n ) THEN
356*
357* Apply the row interchanges to the other blocks.
358*
359 j2 = min( ju-j+1, kv ) - jb
360 j3 = max( 0, ju-j-kv+1 )
361*
362* Use DLASWP to apply the row interchanges to A12, A22, and
363* A32.
364*
365 CALL dlaswp( j2, ab( kv+1-jb, j+jb ), ldab-1, 1, jb,
366 $ ipiv( j ), 1 )
367*
368* Adjust the pivot indices.
369*
370 DO 90 i = j, j + jb - 1
371 ipiv( i ) = ipiv( i ) + j - 1
372 90 CONTINUE
373*
374* Apply the row interchanges to A13, A23, and A33
375* columnwise.
376*
377 k2 = j - 1 + jb + j2
378 DO 110 i = 1, j3
379 jj = k2 + i
380 DO 100 ii = j + i - 1, j + jb - 1
381 ip = ipiv( ii )
382 IF( ip.NE.ii ) THEN
383 temp = ab( kv+1+ii-jj, jj )
384 ab( kv+1+ii-jj, jj ) = ab( kv+1+ip-jj, jj )
385 ab( kv+1+ip-jj, jj ) = temp
386 END IF
387 100 CONTINUE
388 110 CONTINUE
389*
390* Update the relevant part of the trailing submatrix
391*
392 IF( j2.GT.0 ) THEN
393*
394* Update A12
395*
396 CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit',
397 $ jb, j2, one, ab( kv+1, j ), ldab-1,
398 $ ab( kv+1-jb, j+jb ), ldab-1 )
399*
400 IF( i2.GT.0 ) THEN
401*
402* Update A22
403*
404 CALL dgemm( 'No transpose', 'No transpose', i2, j2,
405 $ jb, -one, ab( kv+1+jb, j ), ldab-1,
406 $ ab( kv+1-jb, j+jb ), ldab-1, one,
407 $ ab( kv+1, j+jb ), ldab-1 )
408 END IF
409*
410 IF( i3.GT.0 ) THEN
411*
412* Update A32
413*
414 CALL dgemm( 'No transpose', 'No transpose', i3, j2,
415 $ jb, -one, work31, ldwork,
416 $ ab( kv+1-jb, j+jb ), ldab-1, one,
417 $ ab( kv+kl+1-jb, j+jb ), ldab-1 )
418 END IF
419 END IF
420*
421 IF( j3.GT.0 ) THEN
422*
423* Copy the lower triangle of A13 into the work array
424* WORK13
425*
426 DO 130 jj = 1, j3
427 DO 120 ii = jj, jb
428 work13( ii, jj ) = ab( ii-jj+1, jj+j+kv-1 )
429 120 CONTINUE
430 130 CONTINUE
431*
432* Update A13 in the work array
433*
434 CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit',
435 $ jb, j3, one, ab( kv+1, j ), ldab-1,
436 $ work13, ldwork )
437*
438 IF( i2.GT.0 ) THEN
439*
440* Update A23
441*
442 CALL dgemm( 'No transpose', 'No transpose', i2, j3,
443 $ jb, -one, ab( kv+1+jb, j ), ldab-1,
444 $ work13, ldwork, one, ab( 1+jb, j+kv ),
445 $ ldab-1 )
446 END IF
447*
448 IF( i3.GT.0 ) THEN
449*
450* Update A33
451*
452 CALL dgemm( 'No transpose', 'No transpose', i3, j3,
453 $ jb, -one, work31, ldwork, work13,
454 $ ldwork, one, ab( 1+kl, j+kv ), ldab-1 )
455 END IF
456*
457* Copy the lower triangle of A13 back into place
458*
459 DO 150 jj = 1, j3
460 DO 140 ii = jj, jb
461 ab( ii-jj+1, jj+j+kv-1 ) = work13( ii, jj )
462 140 CONTINUE
463 150 CONTINUE
464 END IF
465 ELSE
466*
467* Adjust the pivot indices.
468*
469 DO 160 i = j, j + jb - 1
470 ipiv( i ) = ipiv( i ) + j - 1
471 160 CONTINUE
472 END IF
473*
474* Partially undo the interchanges in the current block to
475* restore the upper triangular form of A31 and copy the upper
476* triangle of A31 back into place
477*
478 DO 170 jj = j + jb - 1, j, -1
479 jp = ipiv( jj ) - jj + 1
480 IF( jp.NE.1 ) THEN
481*
482* Apply interchange to columns J to JJ-1
483*
484 IF( jp+jj-1.LT.j+kl ) THEN
485*
486* The interchange does not affect A31
487*
488 CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
489 $ ab( kv+jp+jj-j, j ), ldab-1 )
490 ELSE
491*
492* The interchange does affect A31
493*
494 CALL dswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
495 $ work31( jp+jj-j-kl, 1 ), ldwork )
496 END IF
497 END IF
498*
499* Copy the current column of A31 back into place
500*
501 nw = min( i3, jj-j+1 )
502 IF( nw.GT.0 )
503 $ CALL dcopy( nw, work31( 1, jj-j+1 ), 1,
504 $ ab( kv+kl+1-jj+j, jj ), 1 )
505 170 CONTINUE
506 180 CONTINUE
507 END IF
508*
509 RETURN
510*
511* End of DGBTRF
512*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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...
Definition dgbtf2.f:145
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181
Here is the call graph for this function:
Here is the caller graph for this function: