LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dpbtrf.f
Go to the documentation of this file.
1 *> \brief \b DPBTRF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DPBTRF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpbtrf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpbtrf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpbtrf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, KD, LDAB, N
26 * ..
27 * .. Array Arguments ..
28 * DOUBLE PRECISION AB( LDAB, * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> DPBTRF computes the Cholesky factorization of a real symmetric
38 *> positive definite band matrix A.
39 *>
40 *> The factorization has the form
41 *> A = U**T * U, if UPLO = 'U', or
42 *> A = L * L**T, if UPLO = 'L',
43 *> where U is an upper triangular matrix and L is lower triangular.
44 *> \endverbatim
45 *
46 * Arguments:
47 * ==========
48 *
49 *> \param[in] UPLO
50 *> \verbatim
51 *> UPLO is CHARACTER*1
52 *> = 'U': Upper triangle of A is stored;
53 *> = 'L': Lower triangle of A is stored.
54 *> \endverbatim
55 *>
56 *> \param[in] N
57 *> \verbatim
58 *> N is INTEGER
59 *> The order of the matrix A. N >= 0.
60 *> \endverbatim
61 *>
62 *> \param[in] KD
63 *> \verbatim
64 *> KD is INTEGER
65 *> The number of superdiagonals of the matrix A if UPLO = 'U',
66 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
67 *> \endverbatim
68 *>
69 *> \param[in,out] AB
70 *> \verbatim
71 *> AB is DOUBLE PRECISION array, dimension (LDAB,N)
72 *> On entry, the upper or lower triangle of the symmetric band
73 *> matrix A, stored in the first KD+1 rows of the array. The
74 *> j-th column of A is stored in the j-th column of the array AB
75 *> as follows:
76 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
77 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
78 *>
79 *> On exit, if INFO = 0, the triangular factor U or L from the
80 *> Cholesky factorization A = U**T*U or A = L*L**T of the band
81 *> matrix A, in the same storage format as A.
82 *> \endverbatim
83 *>
84 *> \param[in] LDAB
85 *> \verbatim
86 *> LDAB is INTEGER
87 *> The leading dimension of the array AB. LDAB >= KD+1.
88 *> \endverbatim
89 *>
90 *> \param[out] INFO
91 *> \verbatim
92 *> INFO is INTEGER
93 *> = 0: successful exit
94 *> < 0: if INFO = -i, the i-th argument had an illegal value
95 *> > 0: if INFO = i, the leading minor of order i is not
96 *> positive definite, and the factorization could not be
97 *> completed.
98 *> \endverbatim
99 *
100 * Authors:
101 * ========
102 *
103 *> \author Univ. of Tennessee
104 *> \author Univ. of California Berkeley
105 *> \author Univ. of Colorado Denver
106 *> \author NAG Ltd.
107 *
108 *> \date November 2011
109 *
110 *> \ingroup doubleOTHERcomputational
111 *
112 *> \par Further Details:
113 * =====================
114 *>
115 *> \verbatim
116 *>
117 *> The band storage scheme is illustrated by the following example, when
118 *> N = 6, KD = 2, and UPLO = 'U':
119 *>
120 *> On entry: On exit:
121 *>
122 *> * * a13 a24 a35 a46 * * u13 u24 u35 u46
123 *> * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
124 *> a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
125 *>
126 *> Similarly, if UPLO = 'L' the format of A is as follows:
127 *>
128 *> On entry: On exit:
129 *>
130 *> a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66
131 *> a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 *
132 *> a31 a42 a53 a64 * * l31 l42 l53 l64 * *
133 *>
134 *> Array elements marked * are not used by the routine.
135 *> \endverbatim
136 *
137 *> \par Contributors:
138 * ==================
139 *>
140 *> Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
141 *
142 * =====================================================================
143  SUBROUTINE dpbtrf( UPLO, N, KD, AB, LDAB, INFO )
144 *
145 * -- LAPACK computational routine (version 3.4.0) --
146 * -- LAPACK is a software package provided by Univ. of Tennessee, --
147 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148 * November 2011
149 *
150 * .. Scalar Arguments ..
151  CHARACTER uplo
152  INTEGER info, kd, ldab, n
153 * ..
154 * .. Array Arguments ..
155  DOUBLE PRECISION ab( ldab, * )
156 * ..
157 *
158 * =====================================================================
159 *
160 * .. Parameters ..
161  DOUBLE PRECISION one, zero
162  parameter( one = 1.0d+0, zero = 0.0d+0 )
163  INTEGER nbmax, ldwork
164  parameter( nbmax = 32, ldwork = nbmax+1 )
165 * ..
166 * .. Local Scalars ..
167  INTEGER i, i2, i3, ib, ii, j, jj, nb
168 * ..
169 * .. Local Arrays ..
170  DOUBLE PRECISION work( ldwork, nbmax )
171 * ..
172 * .. External Functions ..
173  LOGICAL lsame
174  INTEGER ilaenv
175  EXTERNAL lsame, ilaenv
176 * ..
177 * .. External Subroutines ..
178  EXTERNAL dgemm, dpbtf2, dpotf2, dsyrk, dtrsm, xerbla
179 * ..
180 * .. Intrinsic Functions ..
181  INTRINSIC min
182 * ..
183 * .. Executable Statements ..
184 *
185 * Test the input parameters.
186 *
187  info = 0
188  IF( ( .NOT.lsame( uplo, 'U' ) ) .AND.
189  $ ( .NOT.lsame( uplo, 'L' ) ) ) THEN
190  info = -1
191  ELSE IF( n.LT.0 ) THEN
192  info = -2
193  ELSE IF( kd.LT.0 ) THEN
194  info = -3
195  ELSE IF( ldab.LT.kd+1 ) THEN
196  info = -5
197  END IF
198  IF( info.NE.0 ) THEN
199  CALL xerbla( 'DPBTRF', -info )
200  return
201  END IF
202 *
203 * Quick return if possible
204 *
205  IF( n.EQ.0 )
206  $ return
207 *
208 * Determine the block size for this environment
209 *
210  nb = ilaenv( 1, 'DPBTRF', uplo, n, kd, -1, -1 )
211 *
212 * The block size must not exceed the semi-bandwidth KD, and must not
213 * exceed the limit set by the size of the local array WORK.
214 *
215  nb = min( nb, nbmax )
216 *
217  IF( nb.LE.1 .OR. nb.GT.kd ) THEN
218 *
219 * Use unblocked code
220 *
221  CALL dpbtf2( uplo, n, kd, ab, ldab, info )
222  ELSE
223 *
224 * Use blocked code
225 *
226  IF( lsame( uplo, 'U' ) ) THEN
227 *
228 * Compute the Cholesky factorization of a symmetric band
229 * matrix, given the upper triangle of the matrix in band
230 * storage.
231 *
232 * Zero the upper triangle of the work array.
233 *
234  DO 20 j = 1, nb
235  DO 10 i = 1, j - 1
236  work( i, j ) = zero
237  10 continue
238  20 continue
239 *
240 * Process the band matrix one diagonal block at a time.
241 *
242  DO 70 i = 1, n, nb
243  ib = min( nb, n-i+1 )
244 *
245 * Factorize the diagonal block
246 *
247  CALL dpotf2( uplo, ib, ab( kd+1, i ), ldab-1, ii )
248  IF( ii.NE.0 ) THEN
249  info = i + ii - 1
250  go to 150
251  END IF
252  IF( i+ib.LE.n ) THEN
253 *
254 * Update the relevant part of the trailing submatrix.
255 * If A11 denotes the diagonal block which has just been
256 * factorized, then we need to update the remaining
257 * blocks in the diagram:
258 *
259 * A11 A12 A13
260 * A22 A23
261 * A33
262 *
263 * The numbers of rows and columns in the partitioning
264 * are IB, I2, I3 respectively. The blocks A12, A22 and
265 * A23 are empty if IB = KD. The upper triangle of A13
266 * lies outside the band.
267 *
268  i2 = min( kd-ib, n-i-ib+1 )
269  i3 = min( ib, n-i-kd+1 )
270 *
271  IF( i2.GT.0 ) THEN
272 *
273 * Update A12
274 *
275  CALL dtrsm( 'Left', 'Upper', 'Transpose',
276  $ 'Non-unit', ib, i2, one, ab( kd+1, i ),
277  $ ldab-1, ab( kd+1-ib, i+ib ), ldab-1 )
278 *
279 * Update A22
280 *
281  CALL dsyrk( 'Upper', 'Transpose', i2, ib, -one,
282  $ ab( kd+1-ib, i+ib ), ldab-1, one,
283  $ ab( kd+1, i+ib ), ldab-1 )
284  END IF
285 *
286  IF( i3.GT.0 ) THEN
287 *
288 * Copy the lower triangle of A13 into the work array.
289 *
290  DO 40 jj = 1, i3
291  DO 30 ii = jj, ib
292  work( ii, jj ) = ab( ii-jj+1, jj+i+kd-1 )
293  30 continue
294  40 continue
295 *
296 * Update A13 (in the work array).
297 *
298  CALL dtrsm( 'Left', 'Upper', 'Transpose',
299  $ 'Non-unit', ib, i3, one, ab( kd+1, i ),
300  $ ldab-1, work, ldwork )
301 *
302 * Update A23
303 *
304  IF( i2.GT.0 )
305  $ CALL dgemm( 'Transpose', 'No Transpose', i2, i3,
306  $ ib, -one, ab( kd+1-ib, i+ib ),
307  $ ldab-1, work, ldwork, one,
308  $ ab( 1+ib, i+kd ), ldab-1 )
309 *
310 * Update A33
311 *
312  CALL dsyrk( 'Upper', 'Transpose', i3, ib, -one,
313  $ work, ldwork, one, ab( kd+1, i+kd ),
314  $ ldab-1 )
315 *
316 * Copy the lower triangle of A13 back into place.
317 *
318  DO 60 jj = 1, i3
319  DO 50 ii = jj, ib
320  ab( ii-jj+1, jj+i+kd-1 ) = work( ii, jj )
321  50 continue
322  60 continue
323  END IF
324  END IF
325  70 continue
326  ELSE
327 *
328 * Compute the Cholesky factorization of a symmetric band
329 * matrix, given the lower triangle of the matrix in band
330 * storage.
331 *
332 * Zero the lower triangle of the work array.
333 *
334  DO 90 j = 1, nb
335  DO 80 i = j + 1, nb
336  work( i, j ) = zero
337  80 continue
338  90 continue
339 *
340 * Process the band matrix one diagonal block at a time.
341 *
342  DO 140 i = 1, n, nb
343  ib = min( nb, n-i+1 )
344 *
345 * Factorize the diagonal block
346 *
347  CALL dpotf2( uplo, ib, ab( 1, i ), ldab-1, ii )
348  IF( ii.NE.0 ) THEN
349  info = i + ii - 1
350  go to 150
351  END IF
352  IF( i+ib.LE.n ) THEN
353 *
354 * Update the relevant part of the trailing submatrix.
355 * If A11 denotes the diagonal block which has just been
356 * factorized, then we need to update the remaining
357 * blocks in the diagram:
358 *
359 * A11
360 * A21 A22
361 * A31 A32 A33
362 *
363 * The numbers of rows and columns in the partitioning
364 * are IB, I2, I3 respectively. The blocks A21, A22 and
365 * A32 are empty if IB = KD. The lower triangle of A31
366 * lies outside the band.
367 *
368  i2 = min( kd-ib, n-i-ib+1 )
369  i3 = min( ib, n-i-kd+1 )
370 *
371  IF( i2.GT.0 ) THEN
372 *
373 * Update A21
374 *
375  CALL dtrsm( 'Right', 'Lower', 'Transpose',
376  $ 'Non-unit', i2, ib, one, ab( 1, i ),
377  $ ldab-1, ab( 1+ib, i ), ldab-1 )
378 *
379 * Update A22
380 *
381  CALL dsyrk( 'Lower', 'No Transpose', i2, ib, -one,
382  $ ab( 1+ib, i ), ldab-1, one,
383  $ ab( 1, i+ib ), ldab-1 )
384  END IF
385 *
386  IF( i3.GT.0 ) THEN
387 *
388 * Copy the upper triangle of A31 into the work array.
389 *
390  DO 110 jj = 1, ib
391  DO 100 ii = 1, min( jj, i3 )
392  work( ii, jj ) = ab( kd+1-jj+ii, jj+i-1 )
393  100 continue
394  110 continue
395 *
396 * Update A31 (in the work array).
397 *
398  CALL dtrsm( 'Right', 'Lower', 'Transpose',
399  $ 'Non-unit', i3, ib, one, ab( 1, i ),
400  $ ldab-1, work, ldwork )
401 *
402 * Update A32
403 *
404  IF( i2.GT.0 )
405  $ CALL dgemm( 'No transpose', 'Transpose', i3, i2,
406  $ ib, -one, work, ldwork,
407  $ ab( 1+ib, i ), ldab-1, one,
408  $ ab( 1+kd-ib, i+ib ), ldab-1 )
409 *
410 * Update A33
411 *
412  CALL dsyrk( 'Lower', 'No Transpose', i3, ib, -one,
413  $ work, ldwork, one, ab( 1, i+kd ),
414  $ ldab-1 )
415 *
416 * Copy the upper triangle of A31 back into place.
417 *
418  DO 130 jj = 1, ib
419  DO 120 ii = 1, min( jj, i3 )
420  ab( kd+1-jj+ii, jj+i-1 ) = work( ii, jj )
421  120 continue
422  130 continue
423  END IF
424  END IF
425  140 continue
426  END IF
427  END IF
428  return
429 *
430  150 continue
431  return
432 *
433 * End of DPBTRF
434 *
435  END