LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cpbtrf.f
Go to the documentation of this file.
1*> \brief \b CPBTRF
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CPBTRF + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpbtrf.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpbtrf.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpbtrf.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CPBTRF( UPLO, N, KD, AB, LDAB, INFO )
20*
21* .. Scalar Arguments ..
22* CHARACTER UPLO
23* INTEGER INFO, KD, LDAB, N
24* ..
25* .. Array Arguments ..
26* COMPLEX AB( LDAB, * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> CPBTRF computes the Cholesky factorization of a complex Hermitian
36*> positive definite band matrix A.
37*>
38*> The factorization has the form
39*> A = U**H * U, if UPLO = 'U', or
40*> A = L * L**H, if UPLO = 'L',
41*> where U is an upper triangular matrix and L is lower triangular.
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] UPLO
48*> \verbatim
49*> UPLO is CHARACTER*1
50*> = 'U': Upper triangle of A is stored;
51*> = 'L': Lower triangle of A is stored.
52*> \endverbatim
53*>
54*> \param[in] N
55*> \verbatim
56*> N is INTEGER
57*> The order of the matrix A. N >= 0.
58*> \endverbatim
59*>
60*> \param[in] KD
61*> \verbatim
62*> KD is INTEGER
63*> The number of superdiagonals of the matrix A if UPLO = 'U',
64*> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
65*> \endverbatim
66*>
67*> \param[in,out] AB
68*> \verbatim
69*> AB is COMPLEX array, dimension (LDAB,N)
70*> On entry, the upper or lower triangle of the Hermitian band
71*> matrix A, stored in the first KD+1 rows of the array. The
72*> j-th column of A is stored in the j-th column of the array AB
73*> as follows:
74*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
75*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
76*>
77*> On exit, if INFO = 0, the triangular factor U or L from the
78*> Cholesky factorization A = U**H*U or A = L*L**H of the band
79*> matrix A, in the same storage format as A.
80*> \endverbatim
81*>
82*> \param[in] LDAB
83*> \verbatim
84*> LDAB is INTEGER
85*> The leading dimension of the array AB. LDAB >= KD+1.
86*> \endverbatim
87*>
88*> \param[out] INFO
89*> \verbatim
90*> INFO is INTEGER
91*> = 0: successful exit
92*> < 0: if INFO = -i, the i-th argument had an illegal value
93*> > 0: if INFO = i, the leading principal minor of order i
94*> is not positive, and the factorization could not be
95*> completed.
96*> \endverbatim
97*
98* Authors:
99* ========
100*
101*> \author Univ. of Tennessee
102*> \author Univ. of California Berkeley
103*> \author Univ. of Colorado Denver
104*> \author NAG Ltd.
105*
106*> \ingroup pbtrf
107*
108*> \par Further Details:
109* =====================
110*>
111*> \verbatim
112*>
113*> The band storage scheme is illustrated by the following example, when
114*> N = 6, KD = 2, and UPLO = 'U':
115*>
116*> On entry: On exit:
117*>
118*> * * a13 a24 a35 a46 * * u13 u24 u35 u46
119*> * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
120*> a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
121*>
122*> Similarly, if UPLO = 'L' the format of A is as follows:
123*>
124*> On entry: On exit:
125*>
126*> a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66
127*> a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 *
128*> a31 a42 a53 a64 * * l31 l42 l53 l64 * *
129*>
130*> Array elements marked * are not used by the routine.
131*> \endverbatim
132*
133*> \par Contributors:
134* ==================
135*>
136*> Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
137*
138* =====================================================================
139 SUBROUTINE cpbtrf( UPLO, N, KD, AB, LDAB, INFO )
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 AB( LDAB, * )
151* ..
152*
153* =====================================================================
154*
155* .. Parameters ..
156 REAL ONE, ZERO
157 parameter( one = 1.0e+0, zero = 0.0e+0 )
158 COMPLEX CONE
159 parameter( cone = ( 1.0e+0, 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 COMPLEX WORK( LDWORK, NBMAX )
168* ..
169* .. External Functions ..
170 LOGICAL LSAME
171 INTEGER ILAENV
172 EXTERNAL lsame, ilaenv
173* ..
174* .. External Subroutines ..
175 EXTERNAL cgemm, cherk, cpbtf2, cpotf2, ctrsm,
176 $ xerbla
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( 'CPBTRF', -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, 'CPBTRF', 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 cpbtf2( 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 cpotf2( 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 ctrsm( '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 cherk( '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 ctrsm( '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 cgemm( '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 cherk( '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 cpotf2( 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 ctrsm( '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 cherk( '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 ctrsm( '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 cgemm( '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 cherk( '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 CPBTRF
443*
444 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
Definition cherk.f:173
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:140
subroutine cpbtrf(uplo, n, kd, ab, ldab, info)
CPBTRF
Definition cpbtrf.f:140
subroutine cpotf2(uplo, n, a, lda, info)
CPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblock...
Definition cpotf2.f:107
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180