LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
spbtrf.f
Go to the documentation of this file.
1 *> \brief \b SPBTRF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SPBTRF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spbtrf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spbtrf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spbtrf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SPBTRF( UPLO, N, KD, AB, LDAB, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER UPLO
25 * INTEGER INFO, KD, LDAB, N
26 * ..
27 * .. Array Arguments ..
28 * REAL AB( LDAB, * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> SPBTRF 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 REAL 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 *> \ingroup realOTHERcomputational
109 *
110 *> \par Further Details:
111 * =====================
112 *>
113 *> \verbatim
114 *>
115 *> The band storage scheme is illustrated by the following example, when
116 *> N = 6, KD = 2, and UPLO = 'U':
117 *>
118 *> On entry: On exit:
119 *>
120 *> * * a13 a24 a35 a46 * * u13 u24 u35 u46
121 *> * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
122 *> a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
123 *>
124 *> Similarly, if UPLO = 'L' the format of A is as follows:
125 *>
126 *> On entry: On exit:
127 *>
128 *> a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66
129 *> a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 *
130 *> a31 a42 a53 a64 * * l31 l42 l53 l64 * *
131 *>
132 *> Array elements marked * are not used by the routine.
133 *> \endverbatim
134 *
135 *> \par Contributors:
136 * ==================
137 *>
138 *> Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
139 *
140 * =====================================================================
141  SUBROUTINE spbtrf( UPLO, N, KD, AB, LDAB, INFO )
142 *
143 * -- LAPACK computational routine --
144 * -- LAPACK is a software package provided by Univ. of Tennessee, --
145 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146 *
147 * .. Scalar Arguments ..
148  CHARACTER UPLO
149  INTEGER INFO, KD, LDAB, N
150 * ..
151 * .. Array Arguments ..
152  REAL AB( LDAB, * )
153 * ..
154 *
155 * =====================================================================
156 *
157 * .. Parameters ..
158  REAL ONE, ZERO
159  parameter( one = 1.0e+0, zero = 0.0e+0 )
160  INTEGER NBMAX, LDWORK
161  parameter( nbmax = 32, ldwork = nbmax+1 )
162 * ..
163 * .. Local Scalars ..
164  INTEGER I, I2, I3, IB, II, J, JJ, NB
165 * ..
166 * .. Local Arrays ..
167  REAL WORK( LDWORK, NBMAX )
168 * ..
169 * .. External Functions ..
170  LOGICAL LSAME
171  INTEGER ILAENV
172  EXTERNAL lsame, ilaenv
173 * ..
174 * .. External Subroutines ..
175  EXTERNAL sgemm, spbtf2, spotf2, ssyrk, strsm, xerbla
176 * ..
177 * .. Intrinsic Functions ..
178  INTRINSIC min
179 * ..
180 * .. Executable Statements ..
181 *
182 * Test the input parameters.
183 *
184  info = 0
185  IF( ( .NOT.lsame( uplo, 'U' ) ) .AND.
186  $ ( .NOT.lsame( uplo, 'L' ) ) ) THEN
187  info = -1
188  ELSE IF( n.LT.0 ) THEN
189  info = -2
190  ELSE IF( kd.LT.0 ) THEN
191  info = -3
192  ELSE IF( ldab.LT.kd+1 ) THEN
193  info = -5
194  END IF
195  IF( info.NE.0 ) THEN
196  CALL xerbla( 'SPBTRF', -info )
197  RETURN
198  END IF
199 *
200 * Quick return if possible
201 *
202  IF( n.EQ.0 )
203  $ RETURN
204 *
205 * Determine the block size for this environment
206 *
207  nb = ilaenv( 1, 'SPBTRF', uplo, n, kd, -1, -1 )
208 *
209 * The block size must not exceed the semi-bandwidth KD, and must not
210 * exceed the limit set by the size of the local array WORK.
211 *
212  nb = min( nb, nbmax )
213 *
214  IF( nb.LE.1 .OR. nb.GT.kd ) THEN
215 *
216 * Use unblocked code
217 *
218  CALL spbtf2( uplo, n, kd, ab, ldab, info )
219  ELSE
220 *
221 * Use blocked code
222 *
223  IF( lsame( uplo, 'U' ) ) THEN
224 *
225 * Compute the Cholesky factorization of a symmetric band
226 * matrix, given the upper triangle of the matrix in band
227 * storage.
228 *
229 * Zero the upper triangle of the work array.
230 *
231  DO 20 j = 1, nb
232  DO 10 i = 1, j - 1
233  work( i, j ) = zero
234  10 CONTINUE
235  20 CONTINUE
236 *
237 * Process the band matrix one diagonal block at a time.
238 *
239  DO 70 i = 1, n, nb
240  ib = min( nb, n-i+1 )
241 *
242 * Factorize the diagonal block
243 *
244  CALL spotf2( uplo, ib, ab( kd+1, i ), ldab-1, ii )
245  IF( ii.NE.0 ) THEN
246  info = i + ii - 1
247  GO TO 150
248  END IF
249  IF( i+ib.LE.n ) THEN
250 *
251 * Update the relevant part of the trailing submatrix.
252 * If A11 denotes the diagonal block which has just been
253 * factorized, then we need to update the remaining
254 * blocks in the diagram:
255 *
256 * A11 A12 A13
257 * A22 A23
258 * A33
259 *
260 * The numbers of rows and columns in the partitioning
261 * are IB, I2, I3 respectively. The blocks A12, A22 and
262 * A23 are empty if IB = KD. The upper triangle of A13
263 * lies outside the band.
264 *
265  i2 = min( kd-ib, n-i-ib+1 )
266  i3 = min( ib, n-i-kd+1 )
267 *
268  IF( i2.GT.0 ) THEN
269 *
270 * Update A12
271 *
272  CALL strsm( 'Left', 'Upper', 'Transpose',
273  $ 'Non-unit', ib, i2, one, ab( kd+1, i ),
274  $ ldab-1, ab( kd+1-ib, i+ib ), ldab-1 )
275 *
276 * Update A22
277 *
278  CALL ssyrk( 'Upper', 'Transpose', i2, ib, -one,
279  $ ab( kd+1-ib, i+ib ), ldab-1, one,
280  $ ab( kd+1, i+ib ), ldab-1 )
281  END IF
282 *
283  IF( i3.GT.0 ) THEN
284 *
285 * Copy the lower triangle of A13 into the work array.
286 *
287  DO 40 jj = 1, i3
288  DO 30 ii = jj, ib
289  work( ii, jj ) = ab( ii-jj+1, jj+i+kd-1 )
290  30 CONTINUE
291  40 CONTINUE
292 *
293 * Update A13 (in the work array).
294 *
295  CALL strsm( 'Left', 'Upper', 'Transpose',
296  $ 'Non-unit', ib, i3, one, ab( kd+1, i ),
297  $ ldab-1, work, ldwork )
298 *
299 * Update A23
300 *
301  IF( i2.GT.0 )
302  $ CALL sgemm( 'Transpose', 'No Transpose', i2, i3,
303  $ ib, -one, ab( kd+1-ib, i+ib ),
304  $ ldab-1, work, ldwork, one,
305  $ ab( 1+ib, i+kd ), ldab-1 )
306 *
307 * Update A33
308 *
309  CALL ssyrk( 'Upper', 'Transpose', i3, ib, -one,
310  $ work, ldwork, one, ab( kd+1, i+kd ),
311  $ ldab-1 )
312 *
313 * Copy the lower triangle of A13 back into place.
314 *
315  DO 60 jj = 1, i3
316  DO 50 ii = jj, ib
317  ab( ii-jj+1, jj+i+kd-1 ) = work( ii, jj )
318  50 CONTINUE
319  60 CONTINUE
320  END IF
321  END IF
322  70 CONTINUE
323  ELSE
324 *
325 * Compute the Cholesky factorization of a symmetric band
326 * matrix, given the lower triangle of the matrix in band
327 * storage.
328 *
329 * Zero the lower triangle of the work array.
330 *
331  DO 90 j = 1, nb
332  DO 80 i = j + 1, nb
333  work( i, j ) = zero
334  80 CONTINUE
335  90 CONTINUE
336 *
337 * Process the band matrix one diagonal block at a time.
338 *
339  DO 140 i = 1, n, nb
340  ib = min( nb, n-i+1 )
341 *
342 * Factorize the diagonal block
343 *
344  CALL spotf2( uplo, ib, ab( 1, i ), ldab-1, ii )
345  IF( ii.NE.0 ) THEN
346  info = i + ii - 1
347  GO TO 150
348  END IF
349  IF( i+ib.LE.n ) THEN
350 *
351 * Update the relevant part of the trailing submatrix.
352 * If A11 denotes the diagonal block which has just been
353 * factorized, then we need to update the remaining
354 * blocks in the diagram:
355 *
356 * A11
357 * A21 A22
358 * A31 A32 A33
359 *
360 * The numbers of rows and columns in the partitioning
361 * are IB, I2, I3 respectively. The blocks A21, A22 and
362 * A32 are empty if IB = KD. The lower triangle of A31
363 * lies outside the band.
364 *
365  i2 = min( kd-ib, n-i-ib+1 )
366  i3 = min( ib, n-i-kd+1 )
367 *
368  IF( i2.GT.0 ) THEN
369 *
370 * Update A21
371 *
372  CALL strsm( 'Right', 'Lower', 'Transpose',
373  $ 'Non-unit', i2, ib, one, ab( 1, i ),
374  $ ldab-1, ab( 1+ib, i ), ldab-1 )
375 *
376 * Update A22
377 *
378  CALL ssyrk( 'Lower', 'No Transpose', i2, ib, -one,
379  $ ab( 1+ib, i ), ldab-1, one,
380  $ ab( 1, i+ib ), ldab-1 )
381  END IF
382 *
383  IF( i3.GT.0 ) THEN
384 *
385 * Copy the upper triangle of A31 into the work array.
386 *
387  DO 110 jj = 1, ib
388  DO 100 ii = 1, min( jj, i3 )
389  work( ii, jj ) = ab( kd+1-jj+ii, jj+i-1 )
390  100 CONTINUE
391  110 CONTINUE
392 *
393 * Update A31 (in the work array).
394 *
395  CALL strsm( 'Right', 'Lower', 'Transpose',
396  $ 'Non-unit', i3, ib, one, ab( 1, i ),
397  $ ldab-1, work, ldwork )
398 *
399 * Update A32
400 *
401  IF( i2.GT.0 )
402  $ CALL sgemm( 'No transpose', 'Transpose', i3, i2,
403  $ ib, -one, work, ldwork,
404  $ ab( 1+ib, i ), ldab-1, one,
405  $ ab( 1+kd-ib, i+ib ), ldab-1 )
406 *
407 * Update A33
408 *
409  CALL ssyrk( 'Lower', 'No Transpose', i3, ib, -one,
410  $ work, ldwork, one, ab( 1, i+kd ),
411  $ ldab-1 )
412 *
413 * Copy the upper triangle of A31 back into place.
414 *
415  DO 130 jj = 1, ib
416  DO 120 ii = 1, min( jj, i3 )
417  ab( kd+1-jj+ii, jj+i-1 ) = work( ii, jj )
418  120 CONTINUE
419  130 CONTINUE
420  END IF
421  END IF
422  140 CONTINUE
423  END IF
424  END IF
425  RETURN
426 *
427  150 CONTINUE
428  RETURN
429 *
430 * End of SPBTRF
431 *
432  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine spbtf2(UPLO, N, KD, AB, LDAB, INFO)
SPBTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite band matrix (un...
Definition: spbtf2.f:142
subroutine spbtrf(UPLO, N, KD, AB, LDAB, INFO)
SPBTRF
Definition: spbtrf.f:142
subroutine spotf2(UPLO, N, A, LDA, INFO)
SPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblock...
Definition: spotf2.f:109
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:181
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:169
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187