1: SUBROUTINE DGBTRF( 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: * DGBTRF 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 blocked version of the algorithm, calling Level 3 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 interchanges.
85: *
86: * =====================================================================
87: *
88: * .. Parameters ..
89: DOUBLE PRECISION ONE, ZERO
90: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
91: INTEGER NBMAX, LDWORK
92: PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 )
93: * ..
94: * .. Local Scalars ..
95: INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP,
96: $ JU, K2, KM, KV, NB, NW
97: DOUBLE PRECISION TEMP
98: * ..
99: * .. Local Arrays ..
100: DOUBLE PRECISION WORK13( LDWORK, NBMAX ),
101: $ WORK31( LDWORK, NBMAX )
102: * ..
103: * .. External Functions ..
104: INTEGER IDAMAX, ILAENV
105: EXTERNAL IDAMAX, ILAENV
106: * ..
107: * .. External Subroutines ..
108: EXTERNAL DCOPY, DGBTF2, DGEMM, DGER, DLASWP, DSCAL,
109: $ DSWAP, DTRSM, XERBLA
110: * ..
111: * .. Intrinsic Functions ..
112: INTRINSIC MAX, MIN
113: * ..
114: * .. Executable Statements ..
115: *
116: * KV is the number of superdiagonals in the factor U, allowing for
117: * fill-in
118: *
119: KV = KU + KL
120: *
121: * Test the input parameters.
122: *
123: INFO = 0
124: IF( M.LT.0 ) THEN
125: INFO = -1
126: ELSE IF( N.LT.0 ) THEN
127: INFO = -2
128: ELSE IF( KL.LT.0 ) THEN
129: INFO = -3
130: ELSE IF( KU.LT.0 ) THEN
131: INFO = -4
132: ELSE IF( LDAB.LT.KL+KV+1 ) THEN
133: INFO = -6
134: END IF
135: IF( INFO.NE.0 ) THEN
136: CALL XERBLA( 'DGBTRF', -INFO )
137: RETURN
138: END IF
139: *
140: * Quick return if possible
141: *
142: IF( M.EQ.0 .OR. N.EQ.0 )
143: $ RETURN
144: *
145: * Determine the block size for this environment
146: *
147: NB = ILAENV( 1, 'DGBTRF', ' ', M, N, KL, KU )
148: *
149: * The block size must not exceed the limit set by the size of the
150: * local arrays WORK13 and WORK31.
151: *
152: NB = MIN( NB, NBMAX )
153: *
154: IF( NB.LE.1 .OR. NB.GT.KL ) THEN
155: *
156: * Use unblocked code
157: *
158: CALL DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
159: ELSE
160: *
161: * Use blocked code
162: *
163: * Zero the superdiagonal elements of the work array WORK13
164: *
165: DO 20 J = 1, NB
166: DO 10 I = 1, J - 1
167: WORK13( I, J ) = ZERO
168: 10 CONTINUE
169: 20 CONTINUE
170: *
171: * Zero the subdiagonal elements of the work array WORK31
172: *
173: DO 40 J = 1, NB
174: DO 30 I = J + 1, NB
175: WORK31( I, J ) = ZERO
176: 30 CONTINUE
177: 40 CONTINUE
178: *
179: * Gaussian elimination with partial pivoting
180: *
181: * Set fill-in elements in columns KU+2 to KV to zero
182: *
183: DO 60 J = KU + 2, MIN( KV, N )
184: DO 50 I = KV - J + 2, KL
185: AB( I, J ) = ZERO
186: 50 CONTINUE
187: 60 CONTINUE
188: *
189: * JU is the index of the last column affected by the current
190: * stage of the factorization
191: *
192: JU = 1
193: *
194: DO 180 J = 1, MIN( M, N ), NB
195: JB = MIN( NB, MIN( M, N )-J+1 )
196: *
197: * The active part of the matrix is partitioned
198: *
199: * A11 A12 A13
200: * A21 A22 A23
201: * A31 A32 A33
202: *
203: * Here A11, A21 and A31 denote the current block of JB columns
204: * which is about to be factorized. The number of rows in the
205: * partitioning are JB, I2, I3 respectively, and the numbers
206: * of columns are JB, J2, J3. The superdiagonal elements of A13
207: * and the subdiagonal elements of A31 lie outside the band.
208: *
209: I2 = MIN( KL-JB, M-J-JB+1 )
210: I3 = MIN( JB, M-J-KL+1 )
211: *
212: * J2 and J3 are computed after JU has been updated.
213: *
214: * Factorize the current block of JB columns
215: *
216: DO 80 JJ = J, J + JB - 1
217: *
218: * Set fill-in elements in column JJ+KV to zero
219: *
220: IF( JJ+KV.LE.N ) THEN
221: DO 70 I = 1, KL
222: AB( I, JJ+KV ) = ZERO
223: 70 CONTINUE
224: END IF
225: *
226: * Find pivot and test for singularity. KM is the number of
227: * subdiagonal elements in the current column.
228: *
229: KM = MIN( KL, M-JJ )
230: JP = IDAMAX( KM+1, AB( KV+1, JJ ), 1 )
231: IPIV( JJ ) = JP + JJ - J
232: IF( AB( KV+JP, JJ ).NE.ZERO ) THEN
233: JU = MAX( JU, MIN( JJ+KU+JP-1, N ) )
234: IF( JP.NE.1 ) THEN
235: *
236: * Apply interchange to columns J to J+JB-1
237: *
238: IF( JP+JJ-1.LT.J+KL ) THEN
239: *
240: CALL DSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1,
241: $ AB( KV+JP+JJ-J, J ), LDAB-1 )
242: ELSE
243: *
244: * The interchange affects columns J to JJ-1 of A31
245: * which are stored in the work array WORK31
246: *
247: CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
248: $ WORK31( JP+JJ-J-KL, 1 ), LDWORK )
249: CALL DSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1,
250: $ AB( KV+JP, JJ ), LDAB-1 )
251: END IF
252: END IF
253: *
254: * Compute multipliers
255: *
256: CALL DSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ),
257: $ 1 )
258: *
259: * Update trailing submatrix within the band and within
260: * the current block. JM is the index of the last column
261: * which needs to be updated.
262: *
263: JM = MIN( JU, J+JB-1 )
264: IF( JM.GT.JJ )
265: $ CALL DGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1,
266: $ AB( KV, JJ+1 ), LDAB-1,
267: $ AB( KV+1, JJ+1 ), LDAB-1 )
268: ELSE
269: *
270: * If pivot is zero, set INFO to the index of the pivot
271: * unless a zero pivot has already been found.
272: *
273: IF( INFO.EQ.0 )
274: $ INFO = JJ
275: END IF
276: *
277: * Copy current column of A31 into the work array WORK31
278: *
279: NW = MIN( JJ-J+1, I3 )
280: IF( NW.GT.0 )
281: $ CALL DCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1,
282: $ WORK31( 1, JJ-J+1 ), 1 )
283: 80 CONTINUE
284: IF( J+JB.LE.N ) THEN
285: *
286: * Apply the row interchanges to the other blocks.
287: *
288: J2 = MIN( JU-J+1, KV ) - JB
289: J3 = MAX( 0, JU-J-KV+1 )
290: *
291: * Use DLASWP to apply the row interchanges to A12, A22, and
292: * A32.
293: *
294: CALL DLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB,
295: $ IPIV( J ), 1 )
296: *
297: * Adjust the pivot indices.
298: *
299: DO 90 I = J, J + JB - 1
300: IPIV( I ) = IPIV( I ) + J - 1
301: 90 CONTINUE
302: *
303: * Apply the row interchanges to A13, A23, and A33
304: * columnwise.
305: *
306: K2 = J - 1 + JB + J2
307: DO 110 I = 1, J3
308: JJ = K2 + I
309: DO 100 II = J + I - 1, J + JB - 1
310: IP = IPIV( II )
311: IF( IP.NE.II ) THEN
312: TEMP = AB( KV+1+II-JJ, JJ )
313: AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ )
314: AB( KV+1+IP-JJ, JJ ) = TEMP
315: END IF
316: 100 CONTINUE
317: 110 CONTINUE
318: *
319: * Update the relevant part of the trailing submatrix
320: *
321: IF( J2.GT.0 ) THEN
322: *
323: * Update A12
324: *
325: CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit',
326: $ JB, J2, ONE, AB( KV+1, J ), LDAB-1,
327: $ AB( KV+1-JB, J+JB ), LDAB-1 )
328: *
329: IF( I2.GT.0 ) THEN
330: *
331: * Update A22
332: *
333: CALL DGEMM( 'No transpose', 'No transpose', I2, J2,
334: $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1,
335: $ AB( KV+1-JB, J+JB ), LDAB-1, ONE,
336: $ AB( KV+1, J+JB ), LDAB-1 )
337: END IF
338: *
339: IF( I3.GT.0 ) THEN
340: *
341: * Update A32
342: *
343: CALL DGEMM( 'No transpose', 'No transpose', I3, J2,
344: $ JB, -ONE, WORK31, LDWORK,
345: $ AB( KV+1-JB, J+JB ), LDAB-1, ONE,
346: $ AB( KV+KL+1-JB, J+JB ), LDAB-1 )
347: END IF
348: END IF
349: *
350: IF( J3.GT.0 ) THEN
351: *
352: * Copy the lower triangle of A13 into the work array
353: * WORK13
354: *
355: DO 130 JJ = 1, J3
356: DO 120 II = JJ, JB
357: WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 )
358: 120 CONTINUE
359: 130 CONTINUE
360: *
361: * Update A13 in the work array
362: *
363: CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit',
364: $ JB, J3, ONE, AB( KV+1, J ), LDAB-1,
365: $ WORK13, LDWORK )
366: *
367: IF( I2.GT.0 ) THEN
368: *
369: * Update A23
370: *
371: CALL DGEMM( 'No transpose', 'No transpose', I2, J3,
372: $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1,
373: $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ),
374: $ LDAB-1 )
375: END IF
376: *
377: IF( I3.GT.0 ) THEN
378: *
379: * Update A33
380: *
381: CALL DGEMM( 'No transpose', 'No transpose', I3, J3,
382: $ JB, -ONE, WORK31, LDWORK, WORK13,
383: $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 )
384: END IF
385: *
386: * Copy the lower triangle of A13 back into place
387: *
388: DO 150 JJ = 1, J3
389: DO 140 II = JJ, JB
390: AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ )
391: 140 CONTINUE
392: 150 CONTINUE
393: END IF
394: ELSE
395: *
396: * Adjust the pivot indices.
397: *
398: DO 160 I = J, J + JB - 1
399: IPIV( I ) = IPIV( I ) + J - 1
400: 160 CONTINUE
401: END IF
402: *
403: * Partially undo the interchanges in the current block to
404: * restore the upper triangular form of A31 and copy the upper
405: * triangle of A31 back into place
406: *
407: DO 170 JJ = J + JB - 1, J, -1
408: JP = IPIV( JJ ) - JJ + 1
409: IF( JP.NE.1 ) THEN
410: *
411: * Apply interchange to columns J to JJ-1
412: *
413: IF( JP+JJ-1.LT.J+KL ) THEN
414: *
415: * The interchange does not affect A31
416: *
417: CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
418: $ AB( KV+JP+JJ-J, J ), LDAB-1 )
419: ELSE
420: *
421: * The interchange does affect A31
422: *
423: CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
424: $ WORK31( JP+JJ-J-KL, 1 ), LDWORK )
425: END IF
426: END IF
427: *
428: * Copy the current column of A31 back into place
429: *
430: NW = MIN( I3, JJ-J+1 )
431: IF( NW.GT.0 )
432: $ CALL DCOPY( NW, WORK31( 1, JJ-J+1 ), 1,
433: $ AB( KV+KL+1-JJ+J, JJ ), 1 )
434: 170 CONTINUE
435: 180 CONTINUE
436: END IF
437: *
438: RETURN
439: *
440: * End of DGBTRF
441: *
442: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgbtrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack