LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 subdiaonal 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 complexSYcomputational
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,
196  $ claset, cgbtrf, cgemm, cgetrf,
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 )
396 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
397 c INFO = IINFO+(J+1)*NB
398 c 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 )
591 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
592 c INFO = IINFO+(J+1)*NB
593 c 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 *
658 c CALL CLASWP( J*NB, A( 1, 1 ), LDA,
659 c $ (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)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:180
subroutine cgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
CGBTRF
Definition: cgbtrf.f:144
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 clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
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 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 chetrf_aa_2stage(UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2, WORK, LWORK, INFO)
CHETRF_AA_2STAGE