LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zpbtrf()

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 principal minor of order i
!>                is not positive, and the factorization could not be
!>                completed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 139 of file zpbtrf.f.

140*
141* -- LAPACK computational routine --
142* -- LAPACK is a software package provided by Univ. of Tennessee, --
143* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144*
145* .. Scalar Arguments ..
146 CHARACTER UPLO
147 INTEGER INFO, KD, LDAB, N
148* ..
149* .. Array Arguments ..
150 COMPLEX*16 AB( LDAB, * )
151* ..
152*
153* =====================================================================
154*
155* .. Parameters ..
156 DOUBLE PRECISION ONE, ZERO
157 parameter( one = 1.0d+0, zero = 0.0d+0 )
158 COMPLEX*16 CONE
159 parameter( cone = ( 1.0d+0, 0.0d+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 COMPLEX*16 WORK( LDWORK, NBMAX )
168* ..
169* .. External Functions ..
170 LOGICAL LSAME
171 INTEGER ILAENV
172 EXTERNAL lsame, ilaenv
173* ..
174* .. External Subroutines ..
175 EXTERNAL xerbla, zgemm, zherk, zpbtf2, zpotf2,
176 $ ztrsm
177* ..
178* .. Intrinsic Functions ..
179 INTRINSIC min
180* ..
181* .. Executable Statements ..
182*
183* Test the input parameters.
184*
185 info = 0
186 IF( ( .NOT.lsame( uplo, 'U' ) ) .AND.
187 $ ( .NOT.lsame( uplo, 'L' ) ) ) THEN
188 info = -1
189 ELSE IF( n.LT.0 ) THEN
190 info = -2
191 ELSE IF( kd.LT.0 ) THEN
192 info = -3
193 ELSE IF( ldab.LT.kd+1 ) THEN
194 info = -5
195 END IF
196 IF( info.NE.0 ) THEN
197 CALL xerbla( 'ZPBTRF', -info )
198 RETURN
199 END IF
200*
201* Quick return if possible
202*
203 IF( n.EQ.0 )
204 $ RETURN
205*
206* Determine the block size for this environment
207*
208 nb = ilaenv( 1, 'ZPBTRF', uplo, n, kd, -1, -1 )
209*
210* The block size must not exceed the semi-bandwidth KD, and must not
211* exceed the limit set by the size of the local array WORK.
212*
213 nb = min( nb, nbmax )
214*
215 IF( nb.LE.1 .OR. nb.GT.kd ) THEN
216*
217* Use unblocked code
218*
219 CALL zpbtf2( uplo, n, kd, ab, ldab, info )
220 ELSE
221*
222* Use blocked code
223*
224 IF( lsame( uplo, 'U' ) ) THEN
225*
226* Compute the Cholesky factorization of a Hermitian band
227* matrix, given the upper triangle of the matrix in band
228* storage.
229*
230* Zero the upper triangle of the work array.
231*
232 DO 20 j = 1, nb
233 DO 10 i = 1, j - 1
234 work( i, j ) = zero
235 10 CONTINUE
236 20 CONTINUE
237*
238* Process the band matrix one diagonal block at a time.
239*
240 DO 70 i = 1, n, nb
241 ib = min( nb, n-i+1 )
242*
243* Factorize the diagonal block
244*
245 CALL zpotf2( uplo, ib, ab( kd+1, i ), ldab-1, ii )
246 IF( ii.NE.0 ) THEN
247 info = i + ii - 1
248 GO TO 150
249 END IF
250 IF( i+ib.LE.n ) THEN
251*
252* Update the relevant part of the trailing submatrix.
253* If A11 denotes the diagonal block which has just been
254* factorized, then we need to update the remaining
255* blocks in the diagram:
256*
257* A11 A12 A13
258* A22 A23
259* A33
260*
261* The numbers of rows and columns in the partitioning
262* are IB, I2, I3 respectively. The blocks A12, A22 and
263* A23 are empty if IB = KD. The upper triangle of A13
264* lies outside the band.
265*
266 i2 = min( kd-ib, n-i-ib+1 )
267 i3 = min( ib, n-i-kd+1 )
268*
269 IF( i2.GT.0 ) THEN
270*
271* Update A12
272*
273 CALL ztrsm( 'Left', 'Upper',
274 $ '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 zherk( 'Upper', 'Conjugate transpose', i2,
282 $ ib,
283 $ -one, ab( kd+1-ib, i+ib ), ldab-1, one,
284 $ ab( kd+1, i+ib ), ldab-1 )
285 END IF
286*
287 IF( i3.GT.0 ) THEN
288*
289* Copy the lower triangle of A13 into the work array.
290*
291 DO 40 jj = 1, i3
292 DO 30 ii = jj, ib
293 work( ii, jj ) = ab( ii-jj+1, jj+i+kd-1 )
294 30 CONTINUE
295 40 CONTINUE
296*
297* Update A13 (in the work array).
298*
299 CALL ztrsm( 'Left', 'Upper',
300 $ 'Conjugate transpose',
301 $ 'Non-unit', ib, i3, cone,
302 $ ab( kd+1, i ), ldab-1, work, ldwork )
303*
304* Update A23
305*
306 IF( i2.GT.0 )
307 $ CALL zgemm( 'Conjugate transpose',
308 $ 'No transpose', i2, i3, ib, -cone,
309 $ ab( kd+1-ib, i+ib ), ldab-1, work,
310 $ ldwork, cone, ab( 1+ib, i+kd ),
311 $ ldab-1 )
312*
313* Update A33
314*
315 CALL zherk( 'Upper', 'Conjugate transpose', i3,
316 $ 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,
387 $ -one,
388 $ ab( 1+ib, i ), ldab-1, one,
389 $ ab( 1, i+ib ), ldab-1 )
390 END IF
391*
392 IF( i3.GT.0 ) THEN
393*
394* Copy the upper triangle of A31 into the work array.
395*
396 DO 110 jj = 1, ib
397 DO 100 ii = 1, min( jj, i3 )
398 work( ii, jj ) = ab( kd+1-jj+ii, jj+i-1 )
399 100 CONTINUE
400 110 CONTINUE
401*
402* Update A31 (in the work array).
403*
404 CALL ztrsm( 'Right', 'Lower',
405 $ 'Conjugate transpose', 'Non-unit', i3,
406 $ ib, cone, ab( 1, i ), ldab-1, work,
407 $ ldwork )
408*
409* Update A32
410*
411 IF( i2.GT.0 )
412 $ CALL zgemm( 'No transpose',
413 $ 'Conjugate transpose', i3, i2, ib,
414 $ -cone, work, ldwork, ab( 1+ib, i ),
415 $ ldab-1, cone, ab( 1+kd-ib, i+ib ),
416 $ ldab-1 )
417*
418* Update A33
419*
420 CALL zherk( 'Lower', 'No transpose', i3, ib,
421 $ -one,
422 $ work, ldwork, one, ab( 1, i+kd ),
423 $ ldab-1 )
424*
425* Copy the upper triangle of A31 back into place.
426*
427 DO 130 jj = 1, ib
428 DO 120 ii = 1, min( jj, i3 )
429 ab( kd+1-jj+ii, jj+i-1 ) = work( ii, jj )
430 120 CONTINUE
431 130 CONTINUE
432 END IF
433 END IF
434 140 CONTINUE
435 END IF
436 END IF
437 RETURN
438*
439 150 CONTINUE
440 RETURN
441*
442* End of ZPBTRF
443*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zherk
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
Definition zherk.f:173
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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:140
zpotf2
subroutine zpotf2(uplo, n, a, lda, info)
ZPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblock...
Definition zpotf2.f:107
ztrsm
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180
Here is the call graph for this function:
Here is the caller graph for this function: