LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zhetrf_aa_2stage.f
Go to the documentation of this file.
1 *> \brief \b ZHETRF_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 ZHETRF_AA_2STAGE + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_aa_2stage.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHETRF_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*16 A( LDA, * ), TB( * ), WORK( * )
31 * ..
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZHETRF_AA_2STAGE computes the factorization of a double hermitian matrix A
39 *> using the Aasen's algorithm. The form of the factorization is
40 *>
41 *> A = U**H*T*U or A = L*T*L**H
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*16 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*16 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*16 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 complex16SYcomputational
156 *
157 * =====================================================================
158  SUBROUTINE zhetrf_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*16 A( LDA, * ), TB( * ), WORK( * )
174 * ..
175 *
176 * =====================================================================
177 * .. Parameters ..
178  COMPLEX*16 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*16 PIV
187 * ..
188 * .. External Functions ..
189  LOGICAL LSAME
190  INTEGER ILAENV
191  EXTERNAL lsame, ilaenv
192 * ..
193 * .. External Subroutines ..
194  EXTERNAL xerbla, zcopy, zlacgv, zlacpy,
195  $ zlaset, zgbtrf, zgemm, zgetrf,
196  $ zhegst, zswap, ztrsm
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC dconjg, min, max
200 * ..
201 * .. Executable Statements ..
202 *
203 * Test the input parameters.
204 *
205  info = 0
206  upper = lsame( uplo, 'U' )
207  wquery = ( lwork.EQ.-1 )
208  tquery = ( ltb.EQ.-1 )
209  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
210  info = -1
211  ELSE IF( n.LT.0 ) THEN
212  info = -2
213  ELSE IF( lda.LT.max( 1, n ) ) THEN
214  info = -4
215  ELSE IF ( ltb .LT. 4*n .AND. .NOT.tquery ) THEN
216  info = -6
217  ELSE IF ( lwork .LT. n .AND. .NOT.wquery ) THEN
218  info = -10
219  END IF
220 *
221  IF( info.NE.0 ) THEN
222  CALL xerbla( 'ZHETRF_AA_2STAGE', -info )
223  RETURN
224  END IF
225 *
226 * Answer the query
227 *
228  nb = ilaenv( 1, 'ZHETRF_AA_2STAGE', uplo, n, -1, -1, -1 )
229  IF( info.EQ.0 ) THEN
230  IF( tquery ) THEN
231  tb( 1 ) = (3*nb+1)*n
232  END IF
233  IF( wquery ) THEN
234  work( 1 ) = n*nb
235  END IF
236  END IF
237  IF( tquery .OR. wquery ) THEN
238  RETURN
239  END IF
240 *
241 * Quick return
242 *
243  IF ( n.EQ.0 ) THEN
244  RETURN
245  ENDIF
246 *
247 * Determine the number of the block size
248 *
249  ldtb = ltb/n
250  IF( ldtb .LT. 3*nb+1 ) THEN
251  nb = (ldtb-1)/3
252  END IF
253  IF( lwork .LT. nb*n ) THEN
254  nb = lwork/n
255  END IF
256 *
257 * Determine the number of the block columns
258 *
259  nt = (n+nb-1)/nb
260  td = 2*nb
261  kb = min(nb, n)
262 *
263 * Initialize vectors/matrices
264 *
265  DO j = 1, kb
266  ipiv( j ) = j
267  END DO
268 *
269 * Save NB
270 *
271  tb( 1 ) = nb
272 *
273  IF( upper ) THEN
274 *
275 * .....................................................
276 * Factorize A as U**H*D*U using the upper triangle of A
277 * .....................................................
278 *
279  DO j = 0, nt-1
280 *
281 * Generate Jth column of W and H
282 *
283  kb = min(nb, n-j*nb)
284  DO i = 1, j-1
285  IF( i.EQ.1 ) THEN
286 * H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J)
287  IF( i .EQ. (j-1) ) THEN
288  jb = nb+kb
289  ELSE
290  jb = 2*nb
291  END IF
292  CALL zgemm( 'NoTranspose', 'NoTranspose',
293  $ nb, kb, jb,
294  $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
295  $ a( (i-1)*nb+1, j*nb+1 ), lda,
296  $ zero, work( i*nb+1 ), n )
297  ELSE
298 * 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)
299  IF( i .EQ. (j-1) ) THEN
300  jb = 2*nb+kb
301  ELSE
302  jb = 3*nb
303  END IF
304  CALL zgemm( 'NoTranspose', 'NoTranspose',
305  $ nb, kb, jb,
306  $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
307  $ ldtb-1,
308  $ a( (i-2)*nb+1, j*nb+1 ), lda,
309  $ zero, work( i*nb+1 ), n )
310  END IF
311  END DO
312 *
313 * Compute T(J,J)
314 *
315  CALL zlacpy( 'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
316  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
317  IF( j.GT.1 ) THEN
318 * T(J,J) = U(1:J,J)'*H(1:J)
319  CALL zgemm( 'Conjugate transpose', 'NoTranspose',
320  $ kb, kb, (j-1)*nb,
321  $ -one, a( 1, j*nb+1 ), lda,
322  $ work( nb+1 ), n,
323  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
324 * T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
325  CALL zgemm( 'Conjugate transpose', 'NoTranspose',
326  $ kb, nb, kb,
327  $ one, a( (j-1)*nb+1, j*nb+1 ), lda,
328  $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
329  $ zero, work( 1 ), n )
330  CALL zgemm( 'NoTranspose', 'NoTranspose',
331  $ kb, kb, nb,
332  $ -one, work( 1 ), n,
333  $ a( (j-2)*nb+1, j*nb+1 ), lda,
334  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
335  END IF
336  IF( j.GT.0 ) THEN
337  CALL zhegst( 1, 'Upper', kb,
338  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
339  $ a( (j-1)*nb+1, j*nb+1 ), lda, iinfo )
340  END IF
341 *
342 * Expand T(J,J) into full format
343 *
344  DO i = 1, kb
345  tb( td+1 + (j*nb+i-1)*ldtb )
346  $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
347  DO k = i+1, kb
348  tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
349  $ = dconjg( tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb ) )
350  END DO
351  END DO
352 *
353  IF( j.LT.nt-1 ) THEN
354  IF( j.GT.0 ) THEN
355 *
356 * Compute H(J,J)
357 *
358  IF( j.EQ.1 ) THEN
359  CALL zgemm( 'NoTranspose', 'NoTranspose',
360  $ kb, kb, kb,
361  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
362  $ a( (j-1)*nb+1, j*nb+1 ), lda,
363  $ zero, work( j*nb+1 ), n )
364  ELSE
365  CALL zgemm( 'NoTranspose', 'NoTranspose',
366  $ kb, kb, nb+kb,
367  $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
368  $ ldtb-1,
369  $ a( (j-2)*nb+1, j*nb+1 ), lda,
370  $ zero, work( j*nb+1 ), n )
371  END IF
372 *
373 * Update with the previous column
374 *
375  CALL zgemm( 'Conjugate transpose', 'NoTranspose',
376  $ nb, n-(j+1)*nb, j*nb,
377  $ -one, work( nb+1 ), n,
378  $ a( 1, (j+1)*nb+1 ), lda,
379  $ one, a( j*nb+1, (j+1)*nb+1 ), lda )
380  END IF
381 *
382 * Copy panel to workspace to call ZGETRF
383 *
384  DO k = 1, nb
385  CALL zcopy( n-(j+1)*nb,
386  $ a( j*nb+k, (j+1)*nb+1 ), lda,
387  $ work( 1+(k-1)*n ), 1 )
388  END DO
389 *
390 * Factorize panel
391 *
392  CALL zgetrf( n-(j+1)*nb, nb,
393  $ work, n,
394  $ ipiv( (j+1)*nb+1 ), iinfo )
395 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
396 c INFO = IINFO+(J+1)*NB
397 c END IF
398 *
399 * Copy panel back
400 *
401  DO k = 1, nb
402 *
403 * Copy only L-factor
404 *
405  CALL zcopy( n-k-(j+1)*nb,
406  $ work( k+1+(k-1)*n ), 1,
407  $ a( j*nb+k, (j+1)*nb+k+1 ), lda )
408 *
409 * Transpose U-factor to be copied back into T(J+1, J)
410 *
411  CALL zlacgv( k, work( 1+(k-1)*n ), 1 )
412  END DO
413 *
414 * Compute T(J+1, J), zero out for GEMM update
415 *
416  kb = min(nb, n-(j+1)*nb)
417  CALL zlaset( 'Full', kb, nb, zero, zero,
418  $ tb( td+nb+1 + (j*nb)*ldtb) , ldtb-1 )
419  CALL zlacpy( 'Upper', kb, nb,
420  $ work, n,
421  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
422  IF( j.GT.0 ) THEN
423  CALL ztrsm( 'R', 'U', 'N', 'U', kb, nb, one,
424  $ a( (j-1)*nb+1, j*nb+1 ), lda,
425  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
426  END IF
427 *
428 * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
429 * updates
430 *
431  DO k = 1, nb
432  DO i = 1, kb
433  tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
434  $ = dconjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
435  END DO
436  END DO
437  CALL zlaset( 'Lower', kb, nb, zero, one,
438  $ a( j*nb+1, (j+1)*nb+1), lda )
439 *
440 * Apply pivots to trailing submatrix of A
441 *
442  DO k = 1, kb
443 * > Adjust ipiv
444  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
445 *
446  i1 = (j+1)*nb+k
447  i2 = ipiv( (j+1)*nb+k )
448  IF( i1.NE.i2 ) THEN
449 * > Apply pivots to previous columns of L
450  CALL zswap( k-1, a( (j+1)*nb+1, i1 ), 1,
451  $ a( (j+1)*nb+1, i2 ), 1 )
452 * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
453  IF( i2.GT.(i1+1) ) THEN
454  CALL zswap( i2-i1-1, a( i1, i1+1 ), lda,
455  $ a( i1+1, i2 ), 1 )
456  CALL zlacgv( i2-i1-1, a( i1+1, i2 ), 1 )
457  END IF
458  CALL zlacgv( i2-i1, a( i1, i1+1 ), lda )
459 * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
460  IF( i2.LT.n )
461  $ CALL zswap( n-i2, a( i1, i2+1 ), lda,
462  $ a( i2, i2+1 ), lda )
463 * > Swap A(I1, I1) with A(I2, I2)
464  piv = a( i1, i1 )
465  a( i1, i1 ) = a( i2, i2 )
466  a( i2, i2 ) = piv
467 * > Apply pivots to previous columns of L
468  IF( j.GT.0 ) THEN
469  CALL zswap( j*nb, a( 1, i1 ), 1,
470  $ a( 1, i2 ), 1 )
471  END IF
472  ENDIF
473  END DO
474  END IF
475  END DO
476  ELSE
477 *
478 * .....................................................
479 * Factorize A as L*D*L**H using the lower triangle of A
480 * .....................................................
481 *
482  DO j = 0, nt-1
483 *
484 * Generate Jth column of W and H
485 *
486  kb = min(nb, n-j*nb)
487  DO i = 1, j-1
488  IF( i.EQ.1 ) THEN
489 * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
490  IF( i .EQ. (j-1) ) THEN
491  jb = nb+kb
492  ELSE
493  jb = 2*nb
494  END IF
495  CALL zgemm( 'NoTranspose', 'Conjugate transpose',
496  $ nb, kb, jb,
497  $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
498  $ a( j*nb+1, (i-1)*nb+1 ), lda,
499  $ zero, work( i*nb+1 ), n )
500  ELSE
501 * 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)'
502  IF( i .EQ. (j-1) ) THEN
503  jb = 2*nb+kb
504  ELSE
505  jb = 3*nb
506  END IF
507  CALL zgemm( 'NoTranspose', 'Conjugate transpose',
508  $ nb, kb, jb,
509  $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
510  $ ldtb-1,
511  $ a( j*nb+1, (i-2)*nb+1 ), lda,
512  $ zero, work( i*nb+1 ), n )
513  END IF
514  END DO
515 *
516 * Compute T(J,J)
517 *
518  CALL zlacpy( 'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
519  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
520  IF( j.GT.1 ) THEN
521 * T(J,J) = L(J,1:J)*H(1:J)
522  CALL zgemm( 'NoTranspose', 'NoTranspose',
523  $ kb, kb, (j-1)*nb,
524  $ -one, a( j*nb+1, 1 ), lda,
525  $ work( nb+1 ), n,
526  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
527 * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
528  CALL zgemm( 'NoTranspose', 'NoTranspose',
529  $ kb, nb, kb,
530  $ one, a( j*nb+1, (j-1)*nb+1 ), lda,
531  $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
532  $ zero, work( 1 ), n )
533  CALL zgemm( 'NoTranspose', 'Conjugate transpose',
534  $ kb, kb, nb,
535  $ -one, work( 1 ), n,
536  $ a( j*nb+1, (j-2)*nb+1 ), lda,
537  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
538  END IF
539  IF( j.GT.0 ) THEN
540  CALL zhegst( 1, 'Lower', kb,
541  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
542  $ a( j*nb+1, (j-1)*nb+1 ), lda, iinfo )
543  END IF
544 *
545 * Expand T(J,J) into full format
546 *
547  DO i = 1, kb
548  tb( td+1 + (j*nb+i-1)*ldtb )
549  $ = real( tb( td+1 + (j*nb+i-1)*ldtb ) )
550  DO k = i+1, kb
551  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
552  $ = dconjg( tb( td+(k-i)+1 + (j*nb+i-1)*ldtb ) )
553  END DO
554  END DO
555 *
556  IF( j.LT.nt-1 ) THEN
557  IF( j.GT.0 ) THEN
558 *
559 * Compute H(J,J)
560 *
561  IF( j.EQ.1 ) THEN
562  CALL zgemm( 'NoTranspose', 'Conjugate transpose',
563  $ kb, kb, kb,
564  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
565  $ a( j*nb+1, (j-1)*nb+1 ), lda,
566  $ zero, work( j*nb+1 ), n )
567  ELSE
568  CALL zgemm( 'NoTranspose', 'Conjugate transpose',
569  $ kb, kb, nb+kb,
570  $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
571  $ ldtb-1,
572  $ a( j*nb+1, (j-2)*nb+1 ), lda,
573  $ zero, work( j*nb+1 ), n )
574  END IF
575 *
576 * Update with the previous column
577 *
578  CALL zgemm( 'NoTranspose', 'NoTranspose',
579  $ n-(j+1)*nb, nb, j*nb,
580  $ -one, a( (j+1)*nb+1, 1 ), lda,
581  $ work( nb+1 ), n,
582  $ one, a( (j+1)*nb+1, j*nb+1 ), lda )
583  END IF
584 *
585 * Factorize panel
586 *
587  CALL zgetrf( n-(j+1)*nb, nb,
588  $ a( (j+1)*nb+1, j*nb+1 ), lda,
589  $ ipiv( (j+1)*nb+1 ), iinfo )
590 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
591 c INFO = IINFO+(J+1)*NB
592 c END IF
593 *
594 * Compute T(J+1, J), zero out for GEMM update
595 *
596  kb = min(nb, n-(j+1)*nb)
597  CALL zlaset( 'Full', kb, nb, zero, zero,
598  $ tb( td+nb+1 + (j*nb)*ldtb) , ldtb-1 )
599  CALL zlacpy( 'Upper', kb, nb,
600  $ a( (j+1)*nb+1, j*nb+1 ), lda,
601  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
602  IF( j.GT.0 ) THEN
603  CALL ztrsm( 'R', 'L', 'C', 'U', kb, nb, one,
604  $ a( j*nb+1, (j-1)*nb+1 ), lda,
605  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
606  END IF
607 *
608 * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
609 * updates
610 *
611  DO k = 1, nb
612  DO i = 1, kb
613  tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
614  $ = dconjg( tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb ) )
615  END DO
616  END DO
617  CALL zlaset( 'Upper', kb, nb, zero, one,
618  $ a( (j+1)*nb+1, j*nb+1), lda )
619 *
620 * Apply pivots to trailing submatrix of A
621 *
622  DO k = 1, kb
623 * > Adjust ipiv
624  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
625 *
626  i1 = (j+1)*nb+k
627  i2 = ipiv( (j+1)*nb+k )
628  IF( i1.NE.i2 ) THEN
629 * > Apply pivots to previous columns of L
630  CALL zswap( k-1, a( i1, (j+1)*nb+1 ), lda,
631  $ a( i2, (j+1)*nb+1 ), lda )
632 * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
633  IF( i2.GT.(i1+1) ) THEN
634  CALL zswap( i2-i1-1, a( i1+1, i1 ), 1,
635  $ a( i2, i1+1 ), lda )
636  CALL zlacgv( i2-i1-1, a( i2, i1+1 ), lda )
637  END IF
638  CALL zlacgv( i2-i1, a( i1+1, i1 ), 1 )
639 * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
640  IF( i2.LT.n )
641  $ CALL zswap( n-i2, a( i2+1, i1 ), 1,
642  $ a( i2+1, i2 ), 1 )
643 * > Swap A(I1, I1) with A(I2, I2)
644  piv = a( i1, i1 )
645  a( i1, i1 ) = a( i2, i2 )
646  a( i2, i2 ) = piv
647 * > Apply pivots to previous columns of L
648  IF( j.GT.0 ) THEN
649  CALL zswap( j*nb, a( i1, 1 ), lda,
650  $ a( i2, 1 ), lda )
651  END IF
652  ENDIF
653  END DO
654 *
655 * Apply pivots to previous columns of L
656 *
657 c CALL ZLASWP( J*NB, A( 1, 1 ), LDA,
658 c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
659  END IF
660  END DO
661  END IF
662 *
663 * Factor the band matrix
664  CALL zgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
665 *
666  RETURN
667 *
668 * End of ZHETRF_AA_2STAGE
669 *
670  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:180
subroutine zgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
ZGBTRF
Definition: zgbtrf.f:144
subroutine zhegst(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
ZHEGST
Definition: zhegst.f:128
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:74
subroutine zhetrf_aa_2stage(UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2, WORK, LWORK, INFO)
ZHETRF_AA_2STAGE
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102