LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
chetrf_aa_2stage.f
Go to the documentation of this file.
1*> \brief \b CHETRF_AA_2STAGE
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHETRF_AA_2STAGE + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetrf_aa_2stage.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetrf_aa_2stage.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetrf_aa_2stage.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
22* IPIV2, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER UPLO
26* INTEGER N, LDA, LTB, LWORK, INFO
27* ..
28* .. Array Arguments ..
29* INTEGER IPIV( * ), IPIV2( * )
30* COMPLEX A( LDA, * ), TB( * ), WORK( * )
31* ..
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CHETRF_AA_2STAGE computes the factorization of a real hermitian matrix A
39*> using the Aasen's algorithm. The form of the factorization is
40*>
41*> A = U**T*T*U or A = L*T*L**T
42*>
43*> where U (or L) is a product of permutation and unit upper (lower)
44*> triangular matrices, and T is a hermitian band matrix with the
45*> bandwidth of NB (NB is internally selected and stored in TB( 1 ), and T is
46*> LU factorized with partial pivoting).
47*>
48*> This is the blocked version of the algorithm, calling Level 3 BLAS.
49*> \endverbatim
50*
51* Arguments:
52* ==========
53*
54*> \param[in] UPLO
55*> \verbatim
56*> UPLO is CHARACTER*1
57*> = 'U': Upper triangle of A is stored;
58*> = 'L': Lower triangle of A is stored.
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*> N is INTEGER
64*> The order of the matrix A. N >= 0.
65*> \endverbatim
66*>
67*> \param[in,out] A
68*> \verbatim
69*> A is COMPLEX array, dimension (LDA,N)
70*> On entry, the hermitian matrix A. If UPLO = 'U', the leading
71*> N-by-N upper triangular part of A contains the upper
72*> triangular part of the matrix A, and the strictly lower
73*> triangular part of A is not referenced. If UPLO = 'L', the
74*> leading N-by-N lower triangular part of A contains the lower
75*> triangular part of the matrix A, and the strictly upper
76*> triangular part of A is not referenced.
77*>
78*> On exit, L is stored below (or above) the subdiagonal blocks,
79*> when UPLO is 'L' (or 'U').
80*> \endverbatim
81*>
82*> \param[in] LDA
83*> \verbatim
84*> LDA is INTEGER
85*> The leading dimension of the array A. LDA >= max(1,N).
86*> \endverbatim
87*>
88*> \param[out] TB
89*> \verbatim
90*> TB is COMPLEX array, dimension (LTB)
91*> On exit, details of the LU factorization of the band matrix.
92*> \endverbatim
93*>
94*> \param[in] LTB
95*> \verbatim
96*> LTB is INTEGER
97*> The size of the array TB. LTB >= 4*N, internally
98*> used to select NB such that LTB >= (3*NB+1)*N.
99*>
100*> If LTB = -1, then a workspace query is assumed; the
101*> routine only calculates the optimal size of LTB,
102*> returns this value as the first entry of TB, and
103*> no error message related to LTB is issued by XERBLA.
104*> \endverbatim
105*>
106*> \param[out] IPIV
107*> \verbatim
108*> IPIV is INTEGER array, dimension (N)
109*> On exit, it contains the details of the interchanges, i.e.,
110*> the row and column k of A were interchanged with the
111*> row and column IPIV(k).
112*> \endverbatim
113*>
114*> \param[out] IPIV2
115*> \verbatim
116*> IPIV2 is INTEGER array, dimension (N)
117*> On exit, it contains the details of the interchanges, i.e.,
118*> the row and column k of T were interchanged with the
119*> row and column IPIV(k).
120*> \endverbatim
121*>
122*> \param[out] WORK
123*> \verbatim
124*> WORK is COMPLEX workspace of size LWORK
125*> \endverbatim
126*>
127*> \param[in] LWORK
128*> \verbatim
129*> LWORK is INTEGER
130*> The size of WORK. LWORK >= N, internally used to select NB
131*> such that LWORK >= N*NB.
132*>
133*> If LWORK = -1, then a workspace query is assumed; the
134*> routine only calculates the optimal size of the WORK array,
135*> returns this value as the first entry of the WORK array, and
136*> no error message related to LWORK is issued by XERBLA.
137*> \endverbatim
138*>
139*> \param[out] INFO
140*> \verbatim
141*> INFO is INTEGER
142*> = 0: successful exit
143*> < 0: if INFO = -i, the i-th argument had an illegal value.
144*> > 0: if INFO = i, band LU factorization failed on i-th column
145*> \endverbatim
146*
147* Authors:
148* ========
149*
150*> \author Univ. of Tennessee
151*> \author Univ. of California Berkeley
152*> \author Univ. of Colorado Denver
153*> \author NAG Ltd.
154*
155*> \ingroup hetrf_aa_2stage
156*
157* =====================================================================
158 SUBROUTINE chetrf_aa_2stage( UPLO, N, A, LDA, TB, LTB, IPIV,
159 $ IPIV2, WORK, LWORK, INFO )
160*
161* -- LAPACK computational routine --
162* -- LAPACK is a software package provided by Univ. of Tennessee, --
163* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164*
165 IMPLICIT NONE
166*
167* .. Scalar Arguments ..
168 CHARACTER UPLO
169 INTEGER N, LDA, LTB, LWORK, INFO
170* ..
171* .. Array Arguments ..
172 INTEGER IPIV( * ), IPIV2( * )
173 COMPLEX A( LDA, * ), TB( * ), WORK( * )
174* ..
175*
176* =====================================================================
177* .. Parameters ..
178 COMPLEX ZERO, ONE
179 parameter( zero = ( 0.0e+0, 0.0e+0 ),
180 $ one = ( 1.0e+0, 0.0e+0 ) )
181*
182* .. Local Scalars ..
183 LOGICAL UPPER, TQUERY, WQUERY
184 INTEGER I, J, K, I1, I2, TD
185 INTEGER LDTB, NB, KB, JB, NT, IINFO
186 COMPLEX PIV
187* ..
188* .. External Functions ..
189 LOGICAL LSAME
190 INTEGER ILAENV
191 EXTERNAL lsame, ilaenv
192
193* ..
194* .. External Subroutines ..
195 EXTERNAL xerbla, ccopy, clacgv, clacpy,
197 $ chegst, cswap, ctrsm
198* ..
199* .. Intrinsic Functions ..
200 INTRINSIC conjg, min, max
201* ..
202* .. Executable Statements ..
203*
204* Test the input parameters.
205*
206 info = 0
207 upper = lsame( uplo, 'U' )
208 wquery = ( lwork.EQ.-1 )
209 tquery = ( ltb.EQ.-1 )
210 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
211 info = -1
212 ELSE IF( n.LT.0 ) THEN
213 info = -2
214 ELSE IF( lda.LT.max( 1, n ) ) THEN
215 info = -4
216 ELSE IF ( ltb .LT. 4*n .AND. .NOT.tquery ) THEN
217 info = -6
218 ELSE IF ( lwork .LT. n .AND. .NOT.wquery ) THEN
219 info = -10
220 END IF
221*
222 IF( info.NE.0 ) THEN
223 CALL xerbla( 'CHETRF_AA_2STAGE', -info )
224 RETURN
225 END IF
226*
227* Answer the query
228*
229 nb = ilaenv( 1, 'CHETRF_AA_2STAGE', uplo, n, -1, -1, -1 )
230 IF( info.EQ.0 ) THEN
231 IF( tquery ) THEN
232 tb( 1 ) = (3*nb+1)*n
233 END IF
234 IF( wquery ) THEN
235 work( 1 ) = n*nb
236 END IF
237 END IF
238 IF( tquery .OR. wquery ) THEN
239 RETURN
240 END IF
241*
242* Quick return
243*
244 IF ( n.EQ.0 ) THEN
245 RETURN
246 ENDIF
247*
248* Determine the number of the block size
249*
250 ldtb = ltb/n
251 IF( ldtb .LT. 3*nb+1 ) THEN
252 nb = (ldtb-1)/3
253 END IF
254 IF( lwork .LT. nb*n ) THEN
255 nb = lwork/n
256 END IF
257*
258* Determine the number of the block columns
259*
260 nt = (n+nb-1)/nb
261 td = 2*nb
262 kb = min(nb, n)
263*
264* Initialize vectors/matrices
265*
266 DO j = 1, kb
267 ipiv( j ) = j
268 END DO
269*
270* Save NB
271*
272 tb( 1 ) = nb
273*
274 IF( upper ) THEN
275*
276* .....................................................
277* Factorize A as U**T*D*U using the upper triangle of A
278* .....................................................
279*
280 DO j = 0, nt-1
281*
282* Generate Jth column of W and H
283*
284 kb = min(nb, n-j*nb)
285 DO i = 1, j-1
286 IF( i.EQ.1 ) THEN
287* H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J)
288 IF( i .EQ. (j-1) ) THEN
289 jb = nb+kb
290 ELSE
291 jb = 2*nb
292 END IF
293 CALL cgemm( 'NoTranspose', 'NoTranspose',
294 $ nb, kb, jb,
295 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
296 $ a( (i-1)*nb+1, j*nb+1 ), lda,
297 $ zero, work( i*nb+1 ), n )
298 ELSE
299* H(I,J) = T(I,I-1)*U(I-1,J) + T(I,I)*U(I,J) + T(I,I+1)*U(I+1,J)
300 IF( i .EQ. (j-1) ) THEN
301 jb = 2*nb+kb
302 ELSE
303 jb = 3*nb
304 END IF
305 CALL cgemm( 'NoTranspose', 'NoTranspose',
306 $ nb, kb, jb,
307 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
308 $ ldtb-1,
309 $ a( (i-2)*nb+1, j*nb+1 ), lda,
310 $ zero, work( i*nb+1 ), n )
311 END IF
312 END DO
313*
314* Compute T(J,J)
315*
316 CALL clacpy( 'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
317 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
318 IF( j.GT.1 ) THEN
319* T(J,J) = U(1:J,J)'*H(1:J)
320 CALL cgemm( 'Conjugate transpose', 'NoTranspose',
321 $ kb, kb, (j-1)*nb,
322 $ -one, a( 1, j*nb+1 ), lda,
323 $ work( nb+1 ), n,
324 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
325* T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
326 CALL cgemm( 'Conjugate transpose', 'NoTranspose',
327 $ kb, nb, kb,
328 $ one, a( (j-1)*nb+1, j*nb+1 ), lda,
329 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
330 $ zero, work( 1 ), n )
331 CALL cgemm( 'NoTranspose', 'NoTranspose',
332 $ kb, kb, nb,
333 $ -one, work( 1 ), n,
334 $ a( (j-2)*nb+1, j*nb+1 ), lda,
335 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
336 END IF
337 IF( j.GT.0 ) THEN
338 CALL chegst( 1, 'Upper', kb,
339 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
340 $ a( (j-1)*nb+1, j*nb+1 ), lda, iinfo )
341 END IF
342*
343* Expand T(J,J) into full format
344*
345 DO i = 1, kb
346 tb( td+1 + (j*nb+i-1)*ldtb )
347 $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
348 DO k = i+1, kb
349 tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
350 $ = conjg( tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb ) )
351 END DO
352 END DO
353*
354 IF( j.LT.nt-1 ) THEN
355 IF( j.GT.0 ) THEN
356*
357* Compute H(J,J)
358*
359 IF( j.EQ.1 ) THEN
360 CALL cgemm( 'NoTranspose', 'NoTranspose',
361 $ kb, kb, kb,
362 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
363 $ a( (j-1)*nb+1, j*nb+1 ), lda,
364 $ zero, work( j*nb+1 ), n )
365 ELSE
366 CALL cgemm( 'NoTranspose', 'NoTranspose',
367 $ kb, kb, nb+kb,
368 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
369 $ ldtb-1,
370 $ a( (j-2)*nb+1, j*nb+1 ), lda,
371 $ zero, work( j*nb+1 ), n )
372 END IF
373*
374* Update with the previous column
375*
376 CALL cgemm( 'Conjugate transpose', 'NoTranspose',
377 $ nb, n-(j+1)*nb, j*nb,
378 $ -one, work( nb+1 ), n,
379 $ a( 1, (j+1)*nb+1 ), lda,
380 $ one, a( j*nb+1, (j+1)*nb+1 ), lda )
381 END IF
382*
383* Copy panel to workspace to call CGETRF
384*
385 DO k = 1, nb
386 CALL ccopy( n-(j+1)*nb,
387 $ a( j*nb+k, (j+1)*nb+1 ), lda,
388 $ work( 1+(k-1)*n ), 1 )
389 END DO
390*
391* Factorize panel
392*
393 CALL cgetrf( n-(j+1)*nb, nb,
394 $ work, n,
395 $ ipiv( (j+1)*nb+1 ), iinfo )
396c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
397c INFO = IINFO+(J+1)*NB
398c END IF
399*
400* Copy panel back
401*
402 DO k = 1, nb
403*
404* Copy only L-factor
405*
406 CALL ccopy( n-k-(j+1)*nb,
407 $ work( k+1+(k-1)*n ), 1,
408 $ a( j*nb+k, (j+1)*nb+k+1 ), lda )
409*
410* Transpose U-factor to be copied back into T(J+1, J)
411*
412 CALL clacgv( k, work( 1+(k-1)*n ), 1 )
413 END DO
414*
415* Compute T(J+1, J), zero out for GEMM update
416*
417 kb = min(nb, n-(j+1)*nb)
418 CALL claset( 'Full', kb, nb, zero, zero,
419 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
420 CALL clacpy( 'Upper', kb, nb,
421 $ work, n,
422 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
423 IF( j.GT.0 ) THEN
424 CALL ctrsm( 'R', 'U', 'N', 'U', kb, nb, one,
425 $ a( (j-1)*nb+1, j*nb+1 ), lda,
426 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
427 END IF
428*
429* Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
430* updates
431*
432 DO k = 1, nb
433 DO i = 1, kb
434 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
435 $ = conjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
436 END DO
437 END DO
438 CALL claset( 'Lower', kb, nb, zero, one,
439 $ a( j*nb+1, (j+1)*nb+1), lda )
440*
441* Apply pivots to trailing submatrix of A
442*
443 DO k = 1, kb
444* > Adjust ipiv
445 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
446*
447 i1 = (j+1)*nb+k
448 i2 = ipiv( (j+1)*nb+k )
449 IF( i1.NE.i2 ) THEN
450* > Apply pivots to previous columns of L
451 CALL cswap( k-1, a( (j+1)*nb+1, i1 ), 1,
452 $ a( (j+1)*nb+1, i2 ), 1 )
453* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
454 IF( i2.GT.(i1+1) ) THEN
455 CALL cswap( i2-i1-1, a( i1, i1+1 ), lda,
456 $ a( i1+1, i2 ), 1 )
457 CALL clacgv( i2-i1-1, a( i1+1, i2 ), 1 )
458 END IF
459 CALL clacgv( i2-i1, a( i1, i1+1 ), lda )
460* > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
461 IF( i2.LT.n )
462 $ CALL cswap( n-i2, a( i1, i2+1 ), lda,
463 $ a( i2, i2+1 ), lda )
464* > Swap A(I1, I1) with A(I2, I2)
465 piv = a( i1, i1 )
466 a( i1, i1 ) = a( i2, i2 )
467 a( i2, i2 ) = piv
468* > Apply pivots to previous columns of L
469 IF( j.GT.0 ) THEN
470 CALL cswap( j*nb, a( 1, i1 ), 1,
471 $ a( 1, i2 ), 1 )
472 END IF
473 ENDIF
474 END DO
475 END IF
476 END DO
477 ELSE
478*
479* .....................................................
480* Factorize A as L*D*L**T using the lower triangle of A
481* .....................................................
482*
483 DO j = 0, nt-1
484*
485* Generate Jth column of W and H
486*
487 kb = min(nb, n-j*nb)
488 DO i = 1, j-1
489 IF( i.EQ.1 ) THEN
490* H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
491 IF( i .EQ. (j-1) ) THEN
492 jb = nb+kb
493 ELSE
494 jb = 2*nb
495 END IF
496 CALL cgemm( 'NoTranspose', 'Conjugate transpose',
497 $ nb, kb, jb,
498 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
499 $ a( j*nb+1, (i-1)*nb+1 ), lda,
500 $ zero, work( i*nb+1 ), n )
501 ELSE
502* H(I,J) = T(I,I-1)*L(J,I-1)' + T(I,I)*L(J,I)' + T(I,I+1)*L(J,I+1)'
503 IF( i .EQ. (j-1) ) THEN
504 jb = 2*nb+kb
505 ELSE
506 jb = 3*nb
507 END IF
508 CALL cgemm( 'NoTranspose', 'Conjugate transpose',
509 $ nb, kb, jb,
510 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
511 $ ldtb-1,
512 $ a( j*nb+1, (i-2)*nb+1 ), lda,
513 $ zero, work( i*nb+1 ), n )
514 END IF
515 END DO
516*
517* Compute T(J,J)
518*
519 CALL clacpy( 'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
520 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
521 IF( j.GT.1 ) THEN
522* T(J,J) = L(J,1:J)*H(1:J)
523 CALL cgemm( 'NoTranspose', 'NoTranspose',
524 $ kb, kb, (j-1)*nb,
525 $ -one, a( j*nb+1, 1 ), lda,
526 $ work( nb+1 ), n,
527 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
528* T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
529 CALL cgemm( 'NoTranspose', 'NoTranspose',
530 $ kb, nb, kb,
531 $ one, a( j*nb+1, (j-1)*nb+1 ), lda,
532 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
533 $ zero, work( 1 ), n )
534 CALL cgemm( 'NoTranspose', 'Conjugate transpose',
535 $ kb, kb, nb,
536 $ -one, work( 1 ), n,
537 $ a( j*nb+1, (j-2)*nb+1 ), lda,
538 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
539 END IF
540 IF( j.GT.0 ) THEN
541 CALL chegst( 1, 'Lower', kb,
542 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
543 $ a( j*nb+1, (j-1)*nb+1 ), lda, iinfo )
544 END IF
545*
546* Expand T(J,J) into full format
547*
548 DO i = 1, kb
549 tb( td+1 + (j*nb+i-1)*ldtb )
550 $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
551 DO k = i+1, kb
552 tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
553 $ = conjg( tb( td+(k-i)+1 + (j*nb+i-1)*ldtb ) )
554 END DO
555 END DO
556*
557 IF( j.LT.nt-1 ) THEN
558 IF( j.GT.0 ) THEN
559*
560* Compute H(J,J)
561*
562 IF( j.EQ.1 ) THEN
563 CALL cgemm( 'NoTranspose', 'Conjugate transpose',
564 $ kb, kb, kb,
565 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
566 $ a( j*nb+1, (j-1)*nb+1 ), lda,
567 $ zero, work( j*nb+1 ), n )
568 ELSE
569 CALL cgemm( 'NoTranspose', 'Conjugate transpose',
570 $ kb, kb, nb+kb,
571 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
572 $ ldtb-1,
573 $ a( j*nb+1, (j-2)*nb+1 ), lda,
574 $ zero, work( j*nb+1 ), n )
575 END IF
576*
577* Update with the previous column
578*
579 CALL cgemm( 'NoTranspose', 'NoTranspose',
580 $ n-(j+1)*nb, nb, j*nb,
581 $ -one, a( (j+1)*nb+1, 1 ), lda,
582 $ work( nb+1 ), n,
583 $ one, a( (j+1)*nb+1, j*nb+1 ), lda )
584 END IF
585*
586* Factorize panel
587*
588 CALL cgetrf( n-(j+1)*nb, nb,
589 $ a( (j+1)*nb+1, j*nb+1 ), lda,
590 $ ipiv( (j+1)*nb+1 ), iinfo )
591c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
592c INFO = IINFO+(J+1)*NB
593c END IF
594*
595* Compute T(J+1, J), zero out for GEMM update
596*
597 kb = min(nb, n-(j+1)*nb)
598 CALL claset( 'Full', kb, nb, zero, zero,
599 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
600 CALL clacpy( 'Upper', kb, nb,
601 $ a( (j+1)*nb+1, j*nb+1 ), lda,
602 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
603 IF( j.GT.0 ) THEN
604 CALL ctrsm( 'R', 'L', 'C', 'U', kb, nb, one,
605 $ a( j*nb+1, (j-1)*nb+1 ), lda,
606 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
607 END IF
608*
609* Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
610* updates
611*
612 DO k = 1, nb
613 DO i = 1, kb
614 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
615 $ = conjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
616 END DO
617 END DO
618 CALL claset( 'Upper', kb, nb, zero, one,
619 $ a( (j+1)*nb+1, j*nb+1), lda )
620*
621* Apply pivots to trailing submatrix of A
622*
623 DO k = 1, kb
624* > Adjust ipiv
625 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
626*
627 i1 = (j+1)*nb+k
628 i2 = ipiv( (j+1)*nb+k )
629 IF( i1.NE.i2 ) THEN
630* > Apply pivots to previous columns of L
631 CALL cswap( k-1, a( i1, (j+1)*nb+1 ), lda,
632 $ a( i2, (j+1)*nb+1 ), lda )
633* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
634 IF( i2.GT.(i1+1) ) THEN
635 CALL cswap( i2-i1-1, a( i1+1, i1 ), 1,
636 $ a( i2, i1+1 ), lda )
637 CALL clacgv( i2-i1-1, a( i2, i1+1 ), lda )
638 END IF
639 CALL clacgv( i2-i1, a( i1+1, i1 ), 1 )
640* > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
641 IF( i2.LT.n )
642 $ CALL cswap( n-i2, a( i2+1, i1 ), 1,
643 $ a( i2+1, i2 ), 1 )
644* > Swap A(I1, I1) with A(I2, I2)
645 piv = a( i1, i1 )
646 a( i1, i1 ) = a( i2, i2 )
647 a( i2, i2 ) = piv
648* > Apply pivots to previous columns of L
649 IF( j.GT.0 ) THEN
650 CALL cswap( j*nb, a( i1, 1 ), lda,
651 $ a( i2, 1 ), lda )
652 END IF
653 ENDIF
654 END DO
655*
656* Apply pivots to previous columns of L
657*
658c CALL CLASWP( J*NB, A( 1, 1 ), LDA,
659c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
660 END IF
661 END DO
662 END IF
663*
664* Factor the band matrix
665 CALL cgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
666*
667 RETURN
668*
669* End of CHETRF_AA_2STAGE
670*
671 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgbtrf(m, n, kl, ku, ab, ldab, ipiv, info)
CGBTRF
Definition cgbtrf.f:144
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cgetrf(m, n, a, lda, ipiv, info)
CGETRF
Definition cgetrf.f:108
subroutine chegst(itype, uplo, n, a, lda, b, ldb, info)
CHEGST
Definition chegst.f:128
subroutine chetrf_aa_2stage(uplo, n, a, lda, tb, ltb, ipiv, ipiv2, work, lwork, info)
CHETRF_AA_2STAGE
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180