LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
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 *> \htmlonly
9 *> Download CSYTRF_AA_2STAGE + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrf_aa_2stage.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CSYTRF_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 *> CSYTRF_AA_2STAGE computes the factorization of a complex symmetric matrix A
39 *> using the Aasen's algorithm. The form of the factorization is
40 *>
41 *> A = U*T*U**T 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 complex symmetric 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 *> The size of the array TB. LTB >= 4*N, internally
97 *> used to select NB such that LTB >= (3*NB+1)*N.
98 *>
99 *> If LTB = -1, then a workspace query is assumed; the
100 *> routine only calculates the optimal size of LTB,
101 *> returns this value as the first entry of TB, and
102 *> no error message related to LTB is issued by XERBLA.
103 *> \endverbatim
104 *>
105 *> \param[out] IPIV
106 *> \verbatim
107 *> IPIV is INTEGER array, dimension (N)
108 *> On exit, it contains the details of the interchanges, i.e.,
109 *> the row and column k of A were interchanged with the
110 *> row and column IPIV(k).
111 *> \endverbatim
112 *>
113 *> \param[out] IPIV2
114 *> \verbatim
115 *> IPIV is INTEGER array, dimension (N)
116 *> On exit, it contains the details of the interchanges, i.e.,
117 *> the row and column k of T were interchanged with the
118 *> row and column IPIV(k).
119 *> \endverbatim
120 *>
121 *> \param[out] WORK
122 *> \verbatim
123 *> WORK is COMPLEX workspace of size LWORK
124 *> \endverbatim
125 *>
126 *> \param[in] LWORK
127 *> \verbatim
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 *> \date November 2017
154 *
155 *> \ingroup complexSYcomputational
156 *
157 * =====================================================================
158  SUBROUTINE csytrf_aa_2stage( UPLO, N, A, LDA, TB, LTB, IPIV,
159  $ IPIV2, WORK, LWORK, INFO )
160 *
161 * -- LAPACK computational routine (version 3.8.0) --
162 * -- LAPACK is a software package provided by Univ. of Tennessee, --
163 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164 * November 2017
165 *
166  IMPLICIT NONE
167 *
168 * .. Scalar Arguments ..
169  CHARACTER UPLO
170  INTEGER N, LDA, LTB, LWORK, INFO
171 * ..
172 * .. Array Arguments ..
173  INTEGER IPIV( * ), IPIV2( * )
174  COMPLEX A( lda, * ), TB( * ), WORK( * )
175 * ..
176 *
177 * =====================================================================
178 * .. Parameters ..
179  COMPLEX CZERO, CONE
180  parameter( czero = ( 0.0e+0, 0.0e+0 ),
181  $ cone = ( 1.0e+0, 0.0e+0 ) )
182 *
183 * .. Local Scalars ..
184  LOGICAL UPPER, TQUERY, WQUERY
185  INTEGER I, J, K, I1, I2, TD
186  INTEGER LDTB, NB, KB, JB, NT, IINFO
187  COMPLEX PIV
188 * ..
189 * .. External Functions ..
190  LOGICAL LSAME
191  INTEGER ILAENV
192  EXTERNAL lsame, ilaenv
193 * ..
194 * .. External Subroutines ..
195  EXTERNAL ccopy, cgbtrf, cgemm, cgetrf, clacpy,
196  $ claset, ctrsm, cswap, xerbla
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC 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( 'CSYTRF_AA_2STAGE', -info )
223  RETURN
224  END IF
225 *
226 * Answer the query
227 *
228  nb = ilaenv( 1, 'CSYTRF_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 L*D*L**T 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 cgemm( 'NoTranspose', 'NoTranspose',
293  $ nb, kb, jb,
294  $ cone, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
295  $ a( (i-1)*nb+1, j*nb+1 ), lda,
296  $ czero, 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 cgemm( 'NoTranspose', 'NoTranspose',
305  $ nb, kb, jb,
306  $ cone, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
307  $ ldtb-1,
308  $ a( (i-2)*nb+1, j*nb+1 ), lda,
309  $ czero, work( i*nb+1 ), n )
310  END IF
311  END DO
312 *
313 * Compute T(J,J)
314 *
315  CALL clacpy( '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 cgemm( 'Transpose', 'NoTranspose',
320  $ kb, kb, (j-1)*nb,
321  $ -cone, a( 1, j*nb+1 ), lda,
322  $ work( nb+1 ), n,
323  $ cone, 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 cgemm( 'Transpose', 'NoTranspose',
326  $ kb, nb, kb,
327  $ cone, a( (j-1)*nb+1, j*nb+1 ), lda,
328  $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
329  $ czero, work( 1 ), n )
330  CALL cgemm( 'NoTranspose', 'NoTranspose',
331  $ kb, kb, nb,
332  $ -cone, work( 1 ), n,
333  $ a( (j-2)*nb+1, j*nb+1 ), lda,
334  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
335  END IF
336 *
337 * Expand T(J,J) into full format
338 *
339  DO i = 1, kb
340  DO k = i+1, kb
341  tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
342  $ = tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
343  END DO
344  END DO
345  IF( j.GT.0 ) THEN
346 c CALL CHEGST( 1, 'Upper', KB,
347 c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
348 c $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
349  CALL ctrsm( 'L', 'U', 'T', 'N', kb, kb, cone,
350  $ a( (j-1)*nb+1, j*nb+1 ), lda,
351  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
352  CALL ctrsm( 'R', 'U', 'N', 'N', kb, kb, cone,
353  $ a( (j-1)*nb+1, j*nb+1 ), lda,
354  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
355  END IF
356 *
357  IF( j.LT.nt-1 ) THEN
358  IF( j.GT.0 ) THEN
359 *
360 * Compute H(J,J)
361 *
362  IF( j.EQ.1 ) THEN
363  CALL cgemm( 'NoTranspose', 'NoTranspose',
364  $ kb, kb, kb,
365  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
366  $ a( (j-1)*nb+1, j*nb+1 ), lda,
367  $ czero, work( j*nb+1 ), n )
368  ELSE
369  CALL cgemm( 'NoTranspose', 'NoTranspose',
370  $ kb, kb, nb+kb,
371  $ cone, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
372  $ ldtb-1,
373  $ a( (j-2)*nb+1, j*nb+1 ), lda,
374  $ czero, work( j*nb+1 ), n )
375  END IF
376 *
377 * Update with the previous column
378 *
379  CALL cgemm( 'Transpose', 'NoTranspose',
380  $ nb, n-(j+1)*nb, j*nb,
381  $ -cone, work( nb+1 ), n,
382  $ a( 1, (j+1)*nb+1 ), lda,
383  $ cone, a( j*nb+1, (j+1)*nb+1 ), lda )
384  END IF
385 *
386 * Copy panel to workspace to call CGETRF
387 *
388  DO k = 1, nb
389  CALL ccopy( n-(j+1)*nb,
390  $ a( j*nb+k, (j+1)*nb+1 ), lda,
391  $ work( 1+(k-1)*n ), 1 )
392  END DO
393 *
394 * Factorize panel
395 *
396  CALL cgetrf( n-(j+1)*nb, nb,
397  $ work, n,
398  $ ipiv( (j+1)*nb+1 ), iinfo )
399 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
400 c INFO = IINFO+(J+1)*NB
401 c END IF
402 *
403 * Copy panel back
404 *
405  DO k = 1, nb
406  CALL ccopy( n-(j+1)*nb,
407  $ work( 1+(k-1)*n ), 1,
408  $ a( j*nb+k, (j+1)*nb+1 ), lda )
409  END DO
410 *
411 * Compute T(J+1, J), zero out for GEMM update
412 *
413  kb = min(nb, n-(j+1)*nb)
414  CALL claset( 'Full', kb, nb, czero, czero,
415  $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
416  CALL clacpy( 'Upper', kb, nb,
417  $ work, n,
418  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
419  IF( j.GT.0 ) THEN
420  CALL ctrsm( 'R', 'U', 'N', 'U', kb, nb, cone,
421  $ a( (j-1)*nb+1, j*nb+1 ), lda,
422  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
423  END IF
424 *
425 * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
426 * updates
427 *
428  DO k = 1, nb
429  DO i = 1, kb
430  tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
431  $ = tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
432  END DO
433  END DO
434  CALL claset( 'Lower', kb, nb, czero, cone,
435  $ a( j*nb+1, (j+1)*nb+1), lda )
436 *
437 * Apply pivots to trailing submatrix of A
438 *
439  DO k = 1, kb
440 * > Adjust ipiv
441  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
442 *
443  i1 = (j+1)*nb+k
444  i2 = ipiv( (j+1)*nb+k )
445  IF( i1.NE.i2 ) THEN
446 * > Apply pivots to previous columns of L
447  CALL cswap( k-1, a( (j+1)*nb+1, i1 ), 1,
448  $ a( (j+1)*nb+1, i2 ), 1 )
449 * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
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  CALL cswap( n-i2, a( i1, i2+1 ), lda,
454  $ a( i2, i2+1 ), lda )
455 * > Swap A(I1, I1) with A(I2, I2)
456  piv = a( i1, i1 )
457  a( i1, i1 ) = a( i2, i2 )
458  a( i2, i2 ) = piv
459 * > Apply pivots to previous columns of L
460  IF( j.GT.0 ) THEN
461  CALL cswap( j*nb, a( 1, i1 ), 1,
462  $ a( 1, i2 ), 1 )
463  END IF
464  ENDIF
465  END DO
466  END IF
467  END DO
468  ELSE
469 *
470 * .....................................................
471 * Factorize A as L*D*L**T using the lower triangle of A
472 * .....................................................
473 *
474  DO j = 0, nt-1
475 *
476 * Generate Jth column of W and H
477 *
478  kb = min(nb, n-j*nb)
479  DO i = 1, j-1
480  IF( i.EQ.1 ) THEN
481 * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
482  IF( i .EQ. (j-1) ) THEN
483  jb = nb+kb
484  ELSE
485  jb = 2*nb
486  END IF
487  CALL cgemm( 'NoTranspose', 'Transpose',
488  $ nb, kb, jb,
489  $ cone, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
490  $ a( j*nb+1, (i-1)*nb+1 ), lda,
491  $ czero, work( i*nb+1 ), n )
492  ELSE
493 * 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)'
494  IF( i .EQ. (j-1) ) THEN
495  jb = 2*nb+kb
496  ELSE
497  jb = 3*nb
498  END IF
499  CALL cgemm( 'NoTranspose', 'Transpose',
500  $ nb, kb, jb,
501  $ cone, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
502  $ ldtb-1,
503  $ a( j*nb+1, (i-2)*nb+1 ), lda,
504  $ czero, work( i*nb+1 ), n )
505  END IF
506  END DO
507 *
508 * Compute T(J,J)
509 *
510  CALL clacpy( 'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
511  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
512  IF( j.GT.1 ) THEN
513 * T(J,J) = L(J,1:J)*H(1:J)
514  CALL cgemm( 'NoTranspose', 'NoTranspose',
515  $ kb, kb, (j-1)*nb,
516  $ -cone, a( j*nb+1, 1 ), lda,
517  $ work( nb+1 ), n,
518  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
519 * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
520  CALL cgemm( 'NoTranspose', 'NoTranspose',
521  $ kb, nb, kb,
522  $ cone, a( j*nb+1, (j-1)*nb+1 ), lda,
523  $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
524  $ czero, work( 1 ), n )
525  CALL cgemm( 'NoTranspose', 'Transpose',
526  $ kb, kb, nb,
527  $ -cone, work( 1 ), n,
528  $ a( j*nb+1, (j-2)*nb+1 ), lda,
529  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
530  END IF
531 *
532 * Expand T(J,J) into full format
533 *
534  DO i = 1, kb
535  DO k = i+1, kb
536  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
537  $ = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
538  END DO
539  END DO
540  IF( j.GT.0 ) THEN
541 c CALL CHEGST( 1, 'Lower', KB,
542 c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
543 c $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
544  CALL ctrsm( 'L', 'L', 'N', 'N', kb, kb, cone,
545  $ a( j*nb+1, (j-1)*nb+1 ), lda,
546  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
547  CALL ctrsm( 'R', 'L', 'T', 'N', kb, kb, cone,
548  $ a( j*nb+1, (j-1)*nb+1 ), lda,
549  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
550  END IF
551 *
552 * Symmetrize T(J,J)
553 *
554  DO i = 1, kb
555  DO k = i+1, kb
556  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
557  $ = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
558  END DO
559  END DO
560 *
561  IF( j.LT.nt-1 ) THEN
562  IF( j.GT.0 ) THEN
563 *
564 * Compute H(J,J)
565 *
566  IF( j.EQ.1 ) THEN
567  CALL cgemm( 'NoTranspose', 'Transpose',
568  $ kb, kb, kb,
569  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
570  $ a( j*nb+1, (j-1)*nb+1 ), lda,
571  $ czero, work( j*nb+1 ), n )
572  ELSE
573  CALL cgemm( 'NoTranspose', 'Transpose',
574  $ kb, kb, nb+kb,
575  $ cone, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
576  $ ldtb-1,
577  $ a( j*nb+1, (j-2)*nb+1 ), lda,
578  $ czero, work( j*nb+1 ), n )
579  END IF
580 *
581 * Update with the previous column
582 *
583  CALL cgemm( 'NoTranspose', 'NoTranspose',
584  $ n-(j+1)*nb, nb, j*nb,
585  $ -cone, a( (j+1)*nb+1, 1 ), lda,
586  $ work( nb+1 ), n,
587  $ cone, a( (j+1)*nb+1, j*nb+1 ), lda )
588  END IF
589 *
590 * Factorize panel
591 *
592  CALL cgetrf( n-(j+1)*nb, nb,
593  $ a( (j+1)*nb+1, j*nb+1 ), lda,
594  $ ipiv( (j+1)*nb+1 ), iinfo )
595 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
596 c INFO = IINFO+(J+1)*NB
597 c END IF
598 *
599 * Compute T(J+1, J), zero out for GEMM update
600 *
601  kb = min(nb, n-(j+1)*nb)
602  CALL claset( 'Full', kb, nb, czero, czero,
603  $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
604  CALL clacpy( 'Upper', kb, nb,
605  $ a( (j+1)*nb+1, j*nb+1 ), lda,
606  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
607  IF( j.GT.0 ) THEN
608  CALL ctrsm( 'R', 'L', 'T', 'U', kb, nb, cone,
609  $ a( j*nb+1, (j-1)*nb+1 ), lda,
610  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
611  END IF
612 *
613 * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
614 * updates
615 *
616  DO k = 1, nb
617  DO i = 1, kb
618  tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb ) =
619  $ tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
620  END DO
621  END DO
622  CALL claset( 'Upper', kb, nb, czero, cone,
623  $ a( (j+1)*nb+1, j*nb+1 ), lda )
624 *
625 * Apply pivots to trailing submatrix of A
626 *
627  DO k = 1, kb
628 * > Adjust ipiv
629  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
630 *
631  i1 = (j+1)*nb+k
632  i2 = ipiv( (j+1)*nb+k )
633  IF( i1.NE.i2 ) THEN
634 * > Apply pivots to previous columns of L
635  CALL cswap( k-1, a( i1, (j+1)*nb+1 ), lda,
636  $ a( i2, (j+1)*nb+1 ), lda )
637 * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
638  CALL cswap( i2-i1-1, a( i1+1, i1 ), 1,
639  $ a( i2, i1+1 ), lda )
640 * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
641  CALL cswap( 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 cswap( 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 CLASWP( 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 cgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
665 *
666 * End of CSYTRF_AA_2STAGE
667 *
668  END
subroutine cgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
CGBTRF
Definition: cgbtrf.f:146
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:108
subroutine csytrf_aa_2stage(UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2, WORK, LWORK, INFO)
CSYTRF_AA_2STAGE
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:182
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine cgetrf(M, N, A, LDA, IPIV, INFO)
CGETRF
Definition: cgetrf.f:110
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189