LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ sgbtrf()

 subroutine sgbtrf ( integer M, integer N, integer KL, integer KU, real, dimension( ldab, * ) AB, integer LDAB, integer, dimension( * ) IPIV, integer INFO )

SGBTRF

Purpose:
``` SGBTRF 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 REAL 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.```
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 sgbtrf.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  REAL AB( LDAB, * )
155 * ..
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160  REAL ONE, ZERO
161  parameter( one = 1.0e+0, zero = 0.0e+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  REAL TEMP
169 * ..
170 * .. Local Arrays ..
171  REAL WORK13( LDWORK, NBMAX ),
172  \$ WORK31( LDWORK, NBMAX )
173 * ..
174 * .. External Functions ..
175  INTEGER ILAENV, ISAMAX
176  EXTERNAL ilaenv, isamax
177 * ..
178 * .. External Subroutines ..
179  EXTERNAL scopy, sgbtf2, sgemm, sger, slaswp, sscal,
180  \$ sswap, strsm, 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( 'SGBTRF', -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, 'SGBTRF', ' ', 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 sgbtf2( 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 = isamax( 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 sswap( 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 sswap( jj-j, ab( kv+1+jj-j, j ), ldab-1,
319  \$ work31( jp+jj-j-kl, 1 ), ldwork )
320  CALL sswap( 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 sscal( 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 sger( 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 scopy( 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 SLASWP to apply the row interchanges to A12, A22, and
363 * A32.
364 *
365  CALL slaswp( 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 strsm( '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 sgemm( '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 sgemm( '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 strsm( '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 sgemm( '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 sgemm( '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 sswap( 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 sswap( 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 scopy( 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 SGBTRF
512 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sgbtf2(M, N, KL, KU, AB, LDAB, IPIV, INFO)
SGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algo...
Definition: sgbtf2.f:145
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
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:130
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:181
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
Here is the call graph for this function:
Here is the caller graph for this function: