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