LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
sgbtf2.f
Go to the documentation of this file.
1 *> \brief \b SGBTF2 computes the LU factorization of a general band matrix using the unblocked version of the algorithm.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SGBTF2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgbtf2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgbtf2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgbtf2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER INFO, KL, KU, LDAB, M, N
25 * ..
26 * .. Array Arguments ..
27 * INTEGER IPIV( * )
28 * REAL AB( LDAB, * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> SGBTF2 computes an LU factorization of a real m-by-n band matrix A
38 *> using partial pivoting with row interchanges.
39 *>
40 *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] M
47 *> \verbatim
48 *> M is INTEGER
49 *> The number of rows of the matrix A. M >= 0.
50 *> \endverbatim
51 *>
52 *> \param[in] N
53 *> \verbatim
54 *> N is INTEGER
55 *> The number of columns of the matrix A. N >= 0.
56 *> \endverbatim
57 *>
58 *> \param[in] KL
59 *> \verbatim
60 *> KL is INTEGER
61 *> The number of subdiagonals within the band of A. KL >= 0.
62 *> \endverbatim
63 *>
64 *> \param[in] KU
65 *> \verbatim
66 *> KU is INTEGER
67 *> The number of superdiagonals within the band of A. KU >= 0.
68 *> \endverbatim
69 *>
70 *> \param[in,out] AB
71 *> \verbatim
72 *> AB is REAL array, dimension (LDAB,N)
73 *> On entry, the matrix A in band storage, in rows KL+1 to
74 *> 2*KL+KU+1; rows 1 to KL of the array need not be set.
75 *> The j-th column of A is stored in the j-th column of the
76 *> array AB as follows:
77 *> AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
78 *>
79 *> On exit, details of the factorization: U is stored as an
80 *> upper triangular band matrix with KL+KU superdiagonals in
81 *> rows 1 to KL+KU+1, and the multipliers used during the
82 *> factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
83 *> See below for further details.
84 *> \endverbatim
85 *>
86 *> \param[in] LDAB
87 *> \verbatim
88 *> LDAB is INTEGER
89 *> The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
90 *> \endverbatim
91 *>
92 *> \param[out] IPIV
93 *> \verbatim
94 *> IPIV is INTEGER array, dimension (min(M,N))
95 *> The pivot indices; for 1 <= i <= min(M,N), row i of the
96 *> matrix was interchanged with row IPIV(i).
97 *> \endverbatim
98 *>
99 *> \param[out] INFO
100 *> \verbatim
101 *> INFO is INTEGER
102 *> = 0: successful exit
103 *> < 0: if INFO = -i, the i-th argument had an illegal value
104 *> > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
105 *> has been completed, but the factor U is exactly
106 *> singular, and division by zero will occur if it is used
107 *> to solve a system of equations.
108 *> \endverbatim
109 *
110 * Authors:
111 * ========
112 *
113 *> \author Univ. of Tennessee
114 *> \author Univ. of California Berkeley
115 *> \author Univ. of Colorado Denver
116 *> \author NAG Ltd.
117 *
118 *> \date September 2012
119 *
120 *> \ingroup realGBcomputational
121 *
122 *> \par Further Details:
123 * =====================
124 *>
125 *> \verbatim
126 *>
127 *> The band storage scheme is illustrated by the following example, when
128 *> M = N = 6, KL = 2, KU = 1:
129 *>
130 *> On entry: On exit:
131 *>
132 *> * * * + + + * * * u14 u25 u36
133 *> * * + + + + * * u13 u24 u35 u46
134 *> * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
135 *> a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
136 *> a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
137 *> a31 a42 a53 a64 * * m31 m42 m53 m64 * *
138 *>
139 *> Array elements marked * are not used by the routine; elements marked
140 *> + need not be set on entry, but are required by the routine to store
141 *> elements of U, because of fill-in resulting from the row
142 *> interchanges.
143 *> \endverbatim
144 *>
145 * =====================================================================
146  SUBROUTINE sgbtf2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
147 *
148 * -- LAPACK computational routine (version 3.4.2) --
149 * -- LAPACK is a software package provided by Univ. of Tennessee, --
150 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151 * September 2012
152 *
153 * .. Scalar Arguments ..
154  INTEGER INFO, KL, KU, LDAB, M, N
155 * ..
156 * .. Array Arguments ..
157  INTEGER IPIV( * )
158  REAL AB( ldab, * )
159 * ..
160 *
161 * =====================================================================
162 *
163 * .. Parameters ..
164  REAL ONE, ZERO
165  parameter ( one = 1.0e+0, zero = 0.0e+0 )
166 * ..
167 * .. Local Scalars ..
168  INTEGER I, J, JP, JU, KM, KV
169 * ..
170 * .. External Functions ..
171  INTEGER ISAMAX
172  EXTERNAL isamax
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL sger, sscal, sswap, xerbla
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC max, min
179 * ..
180 * .. Executable Statements ..
181 *
182 * KV is the number of superdiagonals in the factor U, allowing for
183 * fill-in.
184 *
185  kv = ku + kl
186 *
187 * Test the input parameters.
188 *
189  info = 0
190  IF( m.LT.0 ) THEN
191  info = -1
192  ELSE IF( n.LT.0 ) THEN
193  info = -2
194  ELSE IF( kl.LT.0 ) THEN
195  info = -3
196  ELSE IF( ku.LT.0 ) THEN
197  info = -4
198  ELSE IF( ldab.LT.kl+kv+1 ) THEN
199  info = -6
200  END IF
201  IF( info.NE.0 ) THEN
202  CALL xerbla( 'SGBTF2', -info )
203  RETURN
204  END IF
205 *
206 * Quick return if possible
207 *
208  IF( m.EQ.0 .OR. n.EQ.0 )
209  $ RETURN
210 *
211 * Gaussian elimination with partial pivoting
212 *
213 * Set fill-in elements in columns KU+2 to KV to zero.
214 *
215  DO 20 j = ku + 2, min( kv, n )
216  DO 10 i = kv - j + 2, kl
217  ab( i, j ) = zero
218  10 CONTINUE
219  20 CONTINUE
220 *
221 * JU is the index of the last column affected by the current stage
222 * of the factorization.
223 *
224  ju = 1
225 *
226  DO 40 j = 1, min( m, n )
227 *
228 * Set fill-in elements in column J+KV to zero.
229 *
230  IF( j+kv.LE.n ) THEN
231  DO 30 i = 1, kl
232  ab( i, j+kv ) = zero
233  30 CONTINUE
234  END IF
235 *
236 * Find pivot and test for singularity. KM is the number of
237 * subdiagonal elements in the current column.
238 *
239  km = min( kl, m-j )
240  jp = isamax( km+1, ab( kv+1, j ), 1 )
241  ipiv( j ) = jp + j - 1
242  IF( ab( kv+jp, j ).NE.zero ) THEN
243  ju = max( ju, min( j+ku+jp-1, n ) )
244 *
245 * Apply interchange to columns J to JU.
246 *
247  IF( jp.NE.1 )
248  $ CALL sswap( ju-j+1, ab( kv+jp, j ), ldab-1,
249  $ ab( kv+1, j ), ldab-1 )
250 *
251  IF( km.GT.0 ) THEN
252 *
253 * Compute multipliers.
254 *
255  CALL sscal( km, one / ab( kv+1, j ), ab( kv+2, j ), 1 )
256 *
257 * Update trailing submatrix within the band.
258 *
259  IF( ju.GT.j )
260  $ CALL sger( km, ju-j, -one, ab( kv+2, j ), 1,
261  $ ab( kv, j+1 ), ldab-1, ab( kv+1, j+1 ),
262  $ ldab-1 )
263  END IF
264  ELSE
265 *
266 * If pivot is zero, set INFO to the index of the pivot
267 * unless a zero pivot has already been found.
268 *
269  IF( info.EQ.0 )
270  $ info = j
271  END IF
272  40 CONTINUE
273  RETURN
274 *
275 * End of SGBTF2
276 *
277  END
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:132
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:147
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53