LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zpbtrf ( character  UPLO,
integer  N,
integer  KD,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
integer  INFO 
)

ZPBTRF

Download ZPBTRF + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 ZPBTRF 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*16 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
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 144 of file zpbtrf.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: