LAPACK  3.9.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*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 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 *> 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 *> \date November 2017
156 *
157 *> \ingroup complexSYcomputational
158 *
159 * =====================================================================
160  SUBROUTINE csytrf_aa_2stage( UPLO, N, A, LDA, TB, LTB, IPIV,
161  $ IPIV2, WORK, LWORK, INFO )
162 *
163 * -- LAPACK computational routine (version 3.8.0) --
164 * -- LAPACK is a software package provided by Univ. of Tennessee, --
165 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166 * November 2017
167 *
168  IMPLICIT NONE
169 *
170 * .. Scalar Arguments ..
171  CHARACTER UPLO
172  INTEGER N, LDA, LTB, LWORK, INFO
173 * ..
174 * .. Array Arguments ..
175  INTEGER IPIV( * ), IPIV2( * )
176  COMPLEX A( LDA, * ), TB( * ), WORK( * )
177 * ..
178 *
179 * =====================================================================
180 * .. Parameters ..
181  COMPLEX CZERO, CONE
182  parameter( czero = ( 0.0e+0, 0.0e+0 ),
183  $ cone = ( 1.0e+0, 0.0e+0 ) )
184 *
185 * .. Local Scalars ..
186  LOGICAL UPPER, TQUERY, WQUERY
187  INTEGER I, J, K, I1, I2, TD
188  INTEGER LDTB, NB, KB, JB, NT, IINFO
189  COMPLEX PIV
190 * ..
191 * .. External Functions ..
192  LOGICAL LSAME
193  INTEGER ILAENV
194  EXTERNAL lsame, ilaenv
195 * ..
196 * .. External Subroutines ..
197  EXTERNAL ccopy, cgbtrf, cgemm, cgetrf, clacpy,
198  $ claset, ctrsm, cswap, xerbla
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC min, max
202 * ..
203 * .. Executable Statements ..
204 *
205 * Test the input parameters.
206 *
207  info = 0
208  upper = lsame( uplo, 'U' )
209  wquery = ( lwork.EQ.-1 )
210  tquery = ( ltb.EQ.-1 )
211  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
212  info = -1
213  ELSE IF( n.LT.0 ) THEN
214  info = -2
215  ELSE IF( lda.LT.max( 1, n ) ) THEN
216  info = -4
217  ELSE IF ( ltb .LT. 4*n .AND. .NOT.tquery ) THEN
218  info = -6
219  ELSE IF ( lwork .LT. n .AND. .NOT.wquery ) THEN
220  info = -10
221  END IF
222 *
223  IF( info.NE.0 ) THEN
224  CALL xerbla( 'CSYTRF_AA_2STAGE', -info )
225  RETURN
226  END IF
227 *
228 * Answer the query
229 *
230  nb = ilaenv( 1, 'CSYTRF_AA_2STAGE', uplo, n, -1, -1, -1 )
231  IF( info.EQ.0 ) THEN
232  IF( tquery ) THEN
233  tb( 1 ) = (3*nb+1)*n
234  END IF
235  IF( wquery ) THEN
236  work( 1 ) = n*nb
237  END IF
238  END IF
239  IF( tquery .OR. wquery ) THEN
240  RETURN
241  END IF
242 *
243 * Quick return
244 *
245  IF ( n.EQ.0 ) THEN
246  RETURN
247  ENDIF
248 *
249 * Determine the number of the block size
250 *
251  ldtb = ltb/n
252  IF( ldtb .LT. 3*nb+1 ) THEN
253  nb = (ldtb-1)/3
254  END IF
255  IF( lwork .LT. nb*n ) THEN
256  nb = lwork/n
257  END IF
258 *
259 * Determine the number of the block columns
260 *
261  nt = (n+nb-1)/nb
262  td = 2*nb
263  kb = min(nb, n)
264 *
265 * Initialize vectors/matrices
266 *
267  DO j = 1, kb
268  ipiv( j ) = j
269  END DO
270 *
271 * Save NB
272 *
273  tb( 1 ) = nb
274 *
275  IF( upper ) THEN
276 *
277 * .....................................................
278 * Factorize A as U**T*D*U using the upper triangle of A
279 * .....................................................
280 *
281  DO j = 0, nt-1
282 *
283 * Generate Jth column of W and H
284 *
285  kb = min(nb, n-j*nb)
286  DO i = 1, j-1
287  IF( i.EQ.1 ) THEN
288 * H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J)
289  IF( i .EQ. (j-1) ) THEN
290  jb = nb+kb
291  ELSE
292  jb = 2*nb
293  END IF
294  CALL cgemm( 'NoTranspose', 'NoTranspose',
295  $ nb, kb, jb,
296  $ cone, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
297  $ a( (i-1)*nb+1, j*nb+1 ), lda,
298  $ czero, work( i*nb+1 ), n )
299  ELSE
300 * 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)
301  IF( i .EQ. j-1) THEN
302  jb = 2*nb+kb
303  ELSE
304  jb = 3*nb
305  END IF
306  CALL cgemm( 'NoTranspose', 'NoTranspose',
307  $ nb, kb, jb,
308  $ cone, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
309  $ ldtb-1,
310  $ a( (i-2)*nb+1, j*nb+1 ), lda,
311  $ czero, work( i*nb+1 ), n )
312  END IF
313  END DO
314 *
315 * Compute T(J,J)
316 *
317  CALL clacpy( 'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
318  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
319  IF( j.GT.1 ) THEN
320 * T(J,J) = U(1:J,J)'*H(1:J)
321  CALL cgemm( 'Transpose', 'NoTranspose',
322  $ kb, kb, (j-1)*nb,
323  $ -cone, a( 1, j*nb+1 ), lda,
324  $ work( nb+1 ), n,
325  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
326 * T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
327  CALL cgemm( 'Transpose', 'NoTranspose',
328  $ kb, nb, kb,
329  $ cone, a( (j-1)*nb+1, j*nb+1 ), lda,
330  $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
331  $ czero, work( 1 ), n )
332  CALL cgemm( 'NoTranspose', 'NoTranspose',
333  $ kb, kb, nb,
334  $ -cone, work( 1 ), n,
335  $ a( (j-2)*nb+1, j*nb+1 ), lda,
336  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
337  END IF
338 *
339 * Expand T(J,J) into full format
340 *
341  DO i = 1, kb
342  DO k = i+1, kb
343  tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
344  $ = tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
345  END DO
346  END DO
347  IF( j.GT.0 ) THEN
348 c CALL CHEGST( 1, 'Upper', KB,
349 c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
350 c $ A( (J-1)*NB+1, J*NB+1 ), LDA, IINFO )
351  CALL ctrsm( 'L', 'U', 'T', 'N', kb, kb, cone,
352  $ a( (j-1)*nb+1, j*nb+1 ), lda,
353  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
354  CALL ctrsm( 'R', 'U', 'N', 'N', kb, kb, cone,
355  $ a( (j-1)*nb+1, j*nb+1 ), lda,
356  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
357  END IF
358 *
359  IF( j.LT.nt-1 ) THEN
360  IF( j.GT.0 ) THEN
361 *
362 * Compute H(J,J)
363 *
364  IF( j.EQ.1 ) THEN
365  CALL cgemm( 'NoTranspose', 'NoTranspose',
366  $ kb, kb, kb,
367  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
368  $ a( (j-1)*nb+1, j*nb+1 ), lda,
369  $ czero, work( j*nb+1 ), n )
370  ELSE
371  CALL cgemm( 'NoTranspose', 'NoTranspose',
372  $ kb, kb, nb+kb,
373  $ cone, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
374  $ ldtb-1,
375  $ a( (j-2)*nb+1, j*nb+1 ), lda,
376  $ czero, work( j*nb+1 ), n )
377  END IF
378 *
379 * Update with the previous column
380 *
381  CALL cgemm( 'Transpose', 'NoTranspose',
382  $ nb, n-(j+1)*nb, j*nb,
383  $ -cone, work( nb+1 ), n,
384  $ a( 1, (j+1)*nb+1 ), lda,
385  $ cone, a( j*nb+1, (j+1)*nb+1 ), lda )
386  END IF
387 *
388 * Copy panel to workspace to call CGETRF
389 *
390  DO k = 1, nb
391  CALL ccopy( n-(j+1)*nb,
392  $ a( j*nb+k, (j+1)*nb+1 ), lda,
393  $ work( 1+(k-1)*n ), 1 )
394  END DO
395 *
396 * Factorize panel
397 *
398  CALL cgetrf( n-(j+1)*nb, nb,
399  $ work, n,
400  $ ipiv( (j+1)*nb+1 ), iinfo )
401 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
402 c INFO = IINFO+(J+1)*NB
403 c END IF
404 *
405 * Copy panel back
406 *
407  DO k = 1, nb
408  CALL ccopy( n-(j+1)*nb,
409  $ work( 1+(k-1)*n ), 1,
410  $ a( j*nb+k, (j+1)*nb+1 ), lda )
411  END DO
412 *
413 * Compute T(J+1, J), zero out for GEMM update
414 *
415  kb = min(nb, n-(j+1)*nb)
416  CALL claset( 'Full', kb, nb, czero, czero,
417  $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
418  CALL clacpy( 'Upper', kb, nb,
419  $ work, n,
420  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
421  IF( j.GT.0 ) THEN
422  CALL ctrsm( 'R', 'U', 'N', 'U', kb, nb, cone,
423  $ a( (j-1)*nb+1, j*nb+1 ), lda,
424  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
425  END IF
426 *
427 * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
428 * updates
429 *
430  DO k = 1, nb
431  DO i = 1, kb
432  tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
433  $ = tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
434  END DO
435  END DO
436  CALL claset( 'Lower', kb, nb, czero, cone,
437  $ a( j*nb+1, (j+1)*nb+1), lda )
438 *
439 * Apply pivots to trailing submatrix of A
440 *
441  DO k = 1, kb
442 * > Adjust ipiv
443  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
444 *
445  i1 = (j+1)*nb+k
446  i2 = ipiv( (j+1)*nb+k )
447  IF( i1.NE.i2 ) THEN
448 * > Apply pivots to previous columns of L
449  CALL cswap( k-1, a( (j+1)*nb+1, i1 ), 1,
450  $ a( (j+1)*nb+1, i2 ), 1 )
451 * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
452  IF( i2.GT.(i1+1) )
453  $ CALL cswap( i2-i1-1, a( i1, i1+1 ), lda,
454  $ a( i1+1, i2 ), 1 )
455 * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
456  IF( i2.LT.n )
457  $ CALL cswap( n-i2, a( i1, i2+1 ), lda,
458  $ a( i2, i2+1 ), lda )
459 * > Swap A(I1, I1) with A(I2, I2)
460  piv = a( i1, i1 )
461  a( i1, i1 ) = a( i2, i2 )
462  a( i2, i2 ) = piv
463 * > Apply pivots to previous columns of L
464  IF( j.GT.0 ) THEN
465  CALL cswap( j*nb, a( 1, i1 ), 1,
466  $ a( 1, i2 ), 1 )
467  END IF
468  ENDIF
469  END DO
470  END IF
471  END DO
472  ELSE
473 *
474 * .....................................................
475 * Factorize A as L*D*L**T using the lower triangle of A
476 * .....................................................
477 *
478  DO j = 0, nt-1
479 *
480 * Generate Jth column of W and H
481 *
482  kb = min(nb, n-j*nb)
483  DO i = 1, j-1
484  IF( i.EQ.1 ) THEN
485 * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
486  IF( i .EQ. (j-1) ) THEN
487  jb = nb+kb
488  ELSE
489  jb = 2*nb
490  END IF
491  CALL cgemm( 'NoTranspose', 'Transpose',
492  $ nb, kb, jb,
493  $ cone, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
494  $ a( j*nb+1, (i-1)*nb+1 ), lda,
495  $ czero, work( i*nb+1 ), n )
496  ELSE
497 * 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)'
498  IF( i .EQ. (j-1) ) THEN
499  jb = 2*nb+kb
500  ELSE
501  jb = 3*nb
502  END IF
503  CALL cgemm( 'NoTranspose', 'Transpose',
504  $ nb, kb, jb,
505  $ cone, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
506  $ ldtb-1,
507  $ a( j*nb+1, (i-2)*nb+1 ), lda,
508  $ czero, work( i*nb+1 ), n )
509  END IF
510  END DO
511 *
512 * Compute T(J,J)
513 *
514  CALL clacpy( 'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
515  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
516  IF( j.GT.1 ) THEN
517 * T(J,J) = L(J,1:J)*H(1:J)
518  CALL cgemm( 'NoTranspose', 'NoTranspose',
519  $ kb, kb, (j-1)*nb,
520  $ -cone, a( j*nb+1, 1 ), lda,
521  $ work( nb+1 ), n,
522  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
523 * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
524  CALL cgemm( 'NoTranspose', 'NoTranspose',
525  $ kb, nb, kb,
526  $ cone, a( j*nb+1, (j-1)*nb+1 ), lda,
527  $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
528  $ czero, work( 1 ), n )
529  CALL cgemm( 'NoTranspose', 'Transpose',
530  $ kb, kb, nb,
531  $ -cone, work( 1 ), n,
532  $ a( j*nb+1, (j-2)*nb+1 ), lda,
533  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
534  END IF
535 *
536 * Expand T(J,J) into full format
537 *
538  DO i = 1, kb
539  DO k = i+1, kb
540  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
541  $ = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
542  END DO
543  END DO
544  IF( j.GT.0 ) THEN
545 c CALL CHEGST( 1, 'Lower', KB,
546 c $ TB( TD+1 + (J*NB)*LDTB ), LDTB-1,
547 c $ A( J*NB+1, (J-1)*NB+1 ), LDA, IINFO )
548  CALL ctrsm( 'L', 'L', 'N', 'N', kb, kb, cone,
549  $ a( j*nb+1, (j-1)*nb+1 ), lda,
550  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
551  CALL ctrsm( 'R', 'L', 'T', 'N', kb, kb, cone,
552  $ a( j*nb+1, (j-1)*nb+1 ), lda,
553  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
554  END IF
555 *
556 * Symmetrize T(J,J)
557 *
558  DO i = 1, kb
559  DO k = i+1, kb
560  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
561  $ = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
562  END DO
563  END DO
564 *
565  IF( j.LT.nt-1 ) THEN
566  IF( j.GT.0 ) THEN
567 *
568 * Compute H(J,J)
569 *
570  IF( j.EQ.1 ) THEN
571  CALL cgemm( 'NoTranspose', 'Transpose',
572  $ kb, kb, kb,
573  $ cone, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
574  $ a( j*nb+1, (j-1)*nb+1 ), lda,
575  $ czero, work( j*nb+1 ), n )
576  ELSE
577  CALL cgemm( 'NoTranspose', 'Transpose',
578  $ kb, kb, nb+kb,
579  $ cone, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
580  $ ldtb-1,
581  $ a( j*nb+1, (j-2)*nb+1 ), lda,
582  $ czero, work( j*nb+1 ), n )
583  END IF
584 *
585 * Update with the previous column
586 *
587  CALL cgemm( 'NoTranspose', 'NoTranspose',
588  $ n-(j+1)*nb, nb, j*nb,
589  $ -cone, a( (j+1)*nb+1, 1 ), lda,
590  $ work( nb+1 ), n,
591  $ cone, a( (j+1)*nb+1, j*nb+1 ), lda )
592  END IF
593 *
594 * Factorize panel
595 *
596  CALL cgetrf( n-(j+1)*nb, nb,
597  $ a( (j+1)*nb+1, j*nb+1 ), lda,
598  $ ipiv( (j+1)*nb+1 ), iinfo )
599 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
600 c INFO = IINFO+(J+1)*NB
601 c END IF
602 *
603 * Compute T(J+1, J), zero out for GEMM update
604 *
605  kb = min(nb, n-(j+1)*nb)
606  CALL claset( 'Full', kb, nb, czero, czero,
607  $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
608  CALL clacpy( 'Upper', kb, nb,
609  $ a( (j+1)*nb+1, j*nb+1 ), lda,
610  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
611  IF( j.GT.0 ) THEN
612  CALL ctrsm( 'R', 'L', 'T', 'U', kb, nb, cone,
613  $ a( j*nb+1, (j-1)*nb+1 ), lda,
614  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
615  END IF
616 *
617 * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
618 * updates
619 *
620  DO k = 1, nb
621  DO i = 1, kb
622  tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb ) =
623  $ tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
624  END DO
625  END DO
626  CALL claset( 'Upper', kb, nb, czero, cone,
627  $ a( (j+1)*nb+1, j*nb+1 ), lda )
628 *
629 * Apply pivots to trailing submatrix of A
630 *
631  DO k = 1, kb
632 * > Adjust ipiv
633  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
634 *
635  i1 = (j+1)*nb+k
636  i2 = ipiv( (j+1)*nb+k )
637  IF( i1.NE.i2 ) THEN
638 * > Apply pivots to previous columns of L
639  CALL cswap( k-1, a( i1, (j+1)*nb+1 ), lda,
640  $ a( i2, (j+1)*nb+1 ), lda )
641 * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
642  IF( i2.GT.(i1+1) )
643  $ CALL cswap( i2-i1-1, a( i1+1, i1 ), 1,
644  $ a( i2, i1+1 ), lda )
645 * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
646  IF( i2.LT.n )
647  $ CALL cswap( n-i2, a( i2+1, i1 ), 1,
648  $ a( i2+1, i2 ), 1 )
649 * > Swap A(I1, I1) with A(I2, I2)
650  piv = a( i1, i1 )
651  a( i1, i1 ) = a( i2, i2 )
652  a( i2, i2 ) = piv
653 * > Apply pivots to previous columns of L
654  IF( j.GT.0 ) THEN
655  CALL cswap( j*nb, a( i1, 1 ), lda,
656  $ a( i2, 1 ), lda )
657  END IF
658  ENDIF
659  END DO
660 *
661 * Apply pivots to previous columns of L
662 *
663 c CALL CLASWP( J*NB, A( 1, 1 ), LDA,
664 c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
665  END IF
666  END DO
667  END IF
668 *
669 * Factor the band matrix
670  CALL cgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
671 *
672  RETURN
673 *
674 * End of CSYTRF_AA_2STAGE
675 *
676  END
cgetrf
subroutine cgetrf(M, N, A, LDA, IPIV, INFO)
CGETRF
Definition: cgetrf.f:110
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
ctrsm
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:182
clacpy
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
claset
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
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:83
csytrf_aa_2stage
subroutine csytrf_aa_2stage(UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2, WORK, LWORK, INFO)
CSYTRF_AA_2STAGE
Definition: csytrf_aa_2stage.f:162
cgbtrf
subroutine cgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
CGBTRF
Definition: cgbtrf.f:146