LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ cpbtrf()

 subroutine cpbtrf ( character UPLO, integer N, integer KD, complex, dimension( ldab, * ) AB, integer LDAB, integer INFO )

CPBTRF

Purpose:
``` CPBTRF computes the Cholesky factorization of a complex Hermitian
positive definite band matrix A.

The factorization has the form
A = U**H * U,  if UPLO = 'U', or
A = L  * L**H,  if UPLO = 'L',
where U is an upper triangular matrix and L is lower triangular.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] KD ``` KD is INTEGER The number of superdiagonals of the matrix A if UPLO = 'U', or the number of subdiagonals if UPLO = 'L'. KD >= 0.``` [in,out] AB ``` AB is COMPLEX array, dimension (LDAB,N) On entry, the upper or lower triangle of the Hermitian band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). On exit, if INFO = 0, the triangular factor U or L from the Cholesky factorization A = U**H*U or A = L*L**H of the band matrix A, in the same storage format as A.``` [in] LDAB ``` LDAB is INTEGER The leading dimension of the array AB. LDAB >= KD+1.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the leading minor of order i is not positive definite, and the factorization could not be completed.```
Further Details:
```  The band storage scheme is illustrated by the following example, when
N = 6, KD = 2, and UPLO = 'U':

On entry:                       On exit:

*    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
*   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66

Similarly, if UPLO = 'L' the format of A is as follows:

On entry:                       On exit:

a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *

Array elements marked * are not used by the routine.```
Contributors:
Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989

Definition at line 141 of file cpbtrf.f.

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  COMPLEX AB( LDAB, * )
153 * ..
154 *
155 * =====================================================================
156 *
157 * .. Parameters ..
158  REAL ONE, ZERO
159  parameter( one = 1.0e+0, zero = 0.0e+0 )
160  COMPLEX CONE
161  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
162  INTEGER NBMAX, LDWORK
163  parameter( nbmax = 32, ldwork = nbmax+1 )
164 * ..
165 * .. Local Scalars ..
166  INTEGER I, I2, I3, IB, II, J, JJ, NB
167 * ..
168 * .. Local Arrays ..
169  COMPLEX WORK( LDWORK, NBMAX )
170 * ..
171 * .. External Functions ..
172  LOGICAL LSAME
173  INTEGER ILAENV
174  EXTERNAL lsame, ilaenv
175 * ..
176 * .. External Subroutines ..
177  EXTERNAL cgemm, cherk, cpbtf2, cpotf2, ctrsm, xerbla
178 * ..
179 * .. Intrinsic Functions ..
180  INTRINSIC min
181 * ..
182 * .. Executable Statements ..
183 *
184 * Test the input parameters.
185 *
186  info = 0
187  IF( ( .NOT.lsame( uplo, 'U' ) ) .AND.
188  \$ ( .NOT.lsame( uplo, 'L' ) ) ) THEN
189  info = -1
190  ELSE IF( n.LT.0 ) THEN
191  info = -2
192  ELSE IF( kd.LT.0 ) THEN
193  info = -3
194  ELSE IF( ldab.LT.kd+1 ) THEN
195  info = -5
196  END IF
197  IF( info.NE.0 ) THEN
198  CALL xerbla( 'CPBTRF', -info )
199  RETURN
200  END IF
201 *
202 * Quick return if possible
203 *
204  IF( n.EQ.0 )
205  \$ RETURN
206 *
207 * Determine the block size for this environment
208 *
209  nb = ilaenv( 1, 'CPBTRF', uplo, n, kd, -1, -1 )
210 *
211 * The block size must not exceed the semi-bandwidth KD, and must not
212 * exceed the limit set by the size of the local array WORK.
213 *
214  nb = min( nb, nbmax )
215 *
216  IF( nb.LE.1 .OR. nb.GT.kd ) THEN
217 *
218 * Use unblocked code
219 *
220  CALL cpbtf2( uplo, n, kd, ab, ldab, info )
221  ELSE
222 *
223 * Use blocked code
224 *
225  IF( lsame( uplo, 'U' ) ) THEN
226 *
227 * Compute the Cholesky factorization of a Hermitian band
228 * matrix, given the upper triangle of the matrix in band
229 * storage.
230 *
231 * Zero the upper triangle of the work array.
232 *
233  DO 20 j = 1, nb
234  DO 10 i = 1, j - 1
235  work( i, j ) = zero
236  10 CONTINUE
237  20 CONTINUE
238 *
239 * Process the band matrix one diagonal block at a time.
240 *
241  DO 70 i = 1, n, nb
242  ib = min( nb, n-i+1 )
243 *
244 * Factorize the diagonal block
245 *
246  CALL cpotf2( uplo, ib, ab( kd+1, i ), ldab-1, ii )
247  IF( ii.NE.0 ) THEN
248  info = i + ii - 1
249  GO TO 150
250  END IF
251  IF( i+ib.LE.n ) THEN
252 *
253 * Update the relevant part of the trailing submatrix.
254 * If A11 denotes the diagonal block which has just been
255 * factorized, then we need to update the remaining
256 * blocks in the diagram:
257 *
258 * A11 A12 A13
259 * A22 A23
260 * A33
261 *
262 * The numbers of rows and columns in the partitioning
263 * are IB, I2, I3 respectively. The blocks A12, A22 and
264 * A23 are empty if IB = KD. The upper triangle of A13
265 * lies outside the band.
266 *
267  i2 = min( kd-ib, n-i-ib+1 )
268  i3 = min( ib, n-i-kd+1 )
269 *
270  IF( i2.GT.0 ) THEN
271 *
272 * Update A12
273 *
274  CALL ctrsm( 'Left', 'Upper', 'Conjugate transpose',
275  \$ 'Non-unit', ib, i2, cone,
276  \$ ab( kd+1, i ), ldab-1,
277  \$ ab( kd+1-ib, i+ib ), ldab-1 )
278 *
279 * Update A22
280 *
281  CALL cherk( 'Upper', 'Conjugate transpose', i2, ib,
282  \$ -one, 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 ctrsm( 'Left', 'Upper', 'Conjugate transpose',
299  \$ 'Non-unit', ib, i3, cone,
300  \$ ab( kd+1, i ), ldab-1, work, ldwork )
301 *
302 * Update A23
303 *
304  IF( i2.GT.0 )
305  \$ CALL cgemm( 'Conjugate transpose',
306  \$ 'No transpose', i2, i3, ib, -cone,
307  \$ ab( kd+1-ib, i+ib ), ldab-1, work,
308  \$ ldwork, cone, ab( 1+ib, i+kd ),
309  \$ ldab-1 )
310 *
311 * Update A33
312 *
313  CALL cherk( 'Upper', 'Conjugate transpose', i3, ib,
314  \$ -one, work, ldwork, one,
315  \$ ab( kd+1, i+kd ), ldab-1 )
316 *
317 * Copy the lower triangle of A13 back into place.
318 *
319  DO 60 jj = 1, i3
320  DO 50 ii = jj, ib
321  ab( ii-jj+1, jj+i+kd-1 ) = work( ii, jj )
322  50 CONTINUE
323  60 CONTINUE
324  END IF
325  END IF
326  70 CONTINUE
327  ELSE
328 *
329 * Compute the Cholesky factorization of a Hermitian band
330 * matrix, given the lower triangle of the matrix in band
331 * storage.
332 *
333 * Zero the lower triangle of the work array.
334 *
335  DO 90 j = 1, nb
336  DO 80 i = j + 1, nb
337  work( i, j ) = zero
338  80 CONTINUE
339  90 CONTINUE
340 *
341 * Process the band matrix one diagonal block at a time.
342 *
343  DO 140 i = 1, n, nb
344  ib = min( nb, n-i+1 )
345 *
346 * Factorize the diagonal block
347 *
348  CALL cpotf2( uplo, ib, ab( 1, i ), ldab-1, ii )
349  IF( ii.NE.0 ) THEN
350  info = i + ii - 1
351  GO TO 150
352  END IF
353  IF( i+ib.LE.n ) THEN
354 *
355 * Update the relevant part of the trailing submatrix.
356 * If A11 denotes the diagonal block which has just been
357 * factorized, then we need to update the remaining
358 * blocks in the diagram:
359 *
360 * A11
361 * A21 A22
362 * A31 A32 A33
363 *
364 * The numbers of rows and columns in the partitioning
365 * are IB, I2, I3 respectively. The blocks A21, A22 and
366 * A32 are empty if IB = KD. The lower triangle of A31
367 * lies outside the band.
368 *
369  i2 = min( kd-ib, n-i-ib+1 )
370  i3 = min( ib, n-i-kd+1 )
371 *
372  IF( i2.GT.0 ) THEN
373 *
374 * Update A21
375 *
376  CALL ctrsm( 'Right', 'Lower',
377  \$ 'Conjugate transpose', 'Non-unit', i2,
378  \$ ib, cone, ab( 1, i ), ldab-1,
379  \$ ab( 1+ib, i ), ldab-1 )
380 *
381 * Update A22
382 *
383  CALL cherk( 'Lower', 'No transpose', i2, ib, -one,
384  \$ ab( 1+ib, i ), ldab-1, one,
385  \$ ab( 1, i+ib ), ldab-1 )
386  END IF
387 *
388  IF( i3.GT.0 ) THEN
389 *
390 * Copy the upper triangle of A31 into the work array.
391 *
392  DO 110 jj = 1, ib
393  DO 100 ii = 1, min( jj, i3 )
394  work( ii, jj ) = ab( kd+1-jj+ii, jj+i-1 )
395  100 CONTINUE
396  110 CONTINUE
397 *
398 * Update A31 (in the work array).
399 *
400  CALL ctrsm( 'Right', 'Lower',
401  \$ 'Conjugate transpose', 'Non-unit', i3,
402  \$ ib, cone, ab( 1, i ), ldab-1, work,
403  \$ ldwork )
404 *
405 * Update A32
406 *
407  IF( i2.GT.0 )
408  \$ CALL cgemm( 'No transpose',
409  \$ 'Conjugate transpose', i3, i2, ib,
410  \$ -cone, work, ldwork, ab( 1+ib, i ),
411  \$ ldab-1, cone, ab( 1+kd-ib, i+ib ),
412  \$ ldab-1 )
413 *
414 * Update A33
415 *
416  CALL cherk( 'Lower', 'No transpose', i3, ib, -one,
417  \$ work, ldwork, one, ab( 1, i+kd ),
418  \$ ldab-1 )
419 *
420 * Copy the upper triangle of A31 back into place.
421 *
422  DO 130 jj = 1, ib
423  DO 120 ii = 1, min( jj, i3 )
424  ab( kd+1-jj+ii, jj+i-1 ) = work( ii, jj )
425  120 CONTINUE
426  130 CONTINUE
427  END IF
428  END IF
429  140 CONTINUE
430  END IF
431  END IF
432  RETURN
433 *
434  150 CONTINUE
435  RETURN
436 *
437 * End of CPBTRF
438 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:173
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:180
subroutine cpbtf2(UPLO, N, KD, AB, LDAB, INFO)
CPBTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite band matrix (un...
Definition: cpbtf2.f:142
subroutine cpotf2(UPLO, N, A, LDA, INFO)
CPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblock...
Definition: cpotf2.f:109
Here is the call graph for this function:
Here is the caller graph for this function: