Annotation of rpl/lapack/lapack/dgbtf2.f, revision 1.7
1.1 bertrand 1: SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, KL, KU, LDAB, M, N
10: * ..
11: * .. Array Arguments ..
12: INTEGER IPIV( * )
13: DOUBLE PRECISION AB( LDAB, * )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * DGBTF2 computes an LU factorization of a real m-by-n band matrix A
20: * using partial pivoting with row interchanges.
21: *
22: * This is the unblocked version of the algorithm, calling Level 2 BLAS.
23: *
24: * Arguments
25: * =========
26: *
27: * M (input) INTEGER
28: * The number of rows of the matrix A. M >= 0.
29: *
30: * N (input) INTEGER
31: * The number of columns of the matrix A. N >= 0.
32: *
33: * KL (input) INTEGER
34: * The number of subdiagonals within the band of A. KL >= 0.
35: *
36: * KU (input) INTEGER
37: * The number of superdiagonals within the band of A. KU >= 0.
38: *
39: * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
40: * On entry, the matrix A in band storage, in rows KL+1 to
41: * 2*KL+KU+1; rows 1 to KL of the array need not be set.
42: * The j-th column of A is stored in the j-th column of the
43: * array AB as follows:
44: * AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
45: *
46: * On exit, details of the factorization: U is stored as an
47: * upper triangular band matrix with KL+KU superdiagonals in
48: * rows 1 to KL+KU+1, and the multipliers used during the
49: * factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
50: * See below for further details.
51: *
52: * LDAB (input) INTEGER
53: * The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
54: *
55: * IPIV (output) INTEGER array, dimension (min(M,N))
56: * The pivot indices; for 1 <= i <= min(M,N), row i of the
57: * matrix was interchanged with row IPIV(i).
58: *
59: * INFO (output) INTEGER
60: * = 0: successful exit
61: * < 0: if INFO = -i, the i-th argument had an illegal value
62: * > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
63: * has been completed, but the factor U is exactly
64: * singular, and division by zero will occur if it is used
65: * to solve a system of equations.
66: *
67: * Further Details
68: * ===============
69: *
70: * The band storage scheme is illustrated by the following example, when
71: * M = N = 6, KL = 2, KU = 1:
72: *
73: * On entry: On exit:
74: *
75: * * * * + + + * * * u14 u25 u36
76: * * * + + + + * * u13 u24 u35 u46
77: * * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
78: * a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
79: * a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
80: * a31 a42 a53 a64 * * m31 m42 m53 m64 * *
81: *
82: * Array elements marked * are not used by the routine; elements marked
83: * + need not be set on entry, but are required by the routine to store
84: * elements of U, because of fill-in resulting from the row
85: * interchanges.
86: *
87: * =====================================================================
88: *
89: * .. Parameters ..
90: DOUBLE PRECISION ONE, ZERO
91: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
92: * ..
93: * .. Local Scalars ..
94: INTEGER I, J, JP, JU, KM, KV
95: * ..
96: * .. External Functions ..
97: INTEGER IDAMAX
98: EXTERNAL IDAMAX
99: * ..
100: * .. External Subroutines ..
101: EXTERNAL DGER, DSCAL, DSWAP, XERBLA
102: * ..
103: * .. Intrinsic Functions ..
104: INTRINSIC MAX, MIN
105: * ..
106: * .. Executable Statements ..
107: *
108: * KV is the number of superdiagonals in the factor U, allowing for
109: * fill-in.
110: *
111: KV = KU + KL
112: *
113: * Test the input parameters.
114: *
115: INFO = 0
116: IF( M.LT.0 ) THEN
117: INFO = -1
118: ELSE IF( N.LT.0 ) THEN
119: INFO = -2
120: ELSE IF( KL.LT.0 ) THEN
121: INFO = -3
122: ELSE IF( KU.LT.0 ) THEN
123: INFO = -4
124: ELSE IF( LDAB.LT.KL+KV+1 ) THEN
125: INFO = -6
126: END IF
127: IF( INFO.NE.0 ) THEN
128: CALL XERBLA( 'DGBTF2', -INFO )
129: RETURN
130: END IF
131: *
132: * Quick return if possible
133: *
134: IF( M.EQ.0 .OR. N.EQ.0 )
135: $ RETURN
136: *
137: * Gaussian elimination with partial pivoting
138: *
139: * Set fill-in elements in columns KU+2 to KV to zero.
140: *
141: DO 20 J = KU + 2, MIN( KV, N )
142: DO 10 I = KV - J + 2, KL
143: AB( I, J ) = ZERO
144: 10 CONTINUE
145: 20 CONTINUE
146: *
147: * JU is the index of the last column affected by the current stage
148: * of the factorization.
149: *
150: JU = 1
151: *
152: DO 40 J = 1, MIN( M, N )
153: *
154: * Set fill-in elements in column J+KV to zero.
155: *
156: IF( J+KV.LE.N ) THEN
157: DO 30 I = 1, KL
158: AB( I, J+KV ) = ZERO
159: 30 CONTINUE
160: END IF
161: *
162: * Find pivot and test for singularity. KM is the number of
163: * subdiagonal elements in the current column.
164: *
165: KM = MIN( KL, M-J )
166: JP = IDAMAX( KM+1, AB( KV+1, J ), 1 )
167: IPIV( J ) = JP + J - 1
168: IF( AB( KV+JP, J ).NE.ZERO ) THEN
169: JU = MAX( JU, MIN( J+KU+JP-1, N ) )
170: *
171: * Apply interchange to columns J to JU.
172: *
173: IF( JP.NE.1 )
174: $ CALL DSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
175: $ AB( KV+1, J ), LDAB-1 )
176: *
177: IF( KM.GT.0 ) THEN
178: *
179: * Compute multipliers.
180: *
181: CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
182: *
183: * Update trailing submatrix within the band.
184: *
185: IF( JU.GT.J )
186: $ CALL DGER( KM, JU-J, -ONE, AB( KV+2, J ), 1,
187: $ AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ),
188: $ LDAB-1 )
189: END IF
190: ELSE
191: *
192: * If pivot is zero, set INFO to the index of the pivot
193: * unless a zero pivot has already been found.
194: *
195: IF( INFO.EQ.0 )
196: $ INFO = J
197: END IF
198: 40 CONTINUE
199: RETURN
200: *
201: * End of DGBTF2
202: *
203: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgbtf2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack