LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
dsytrf_aa_2stage.f
Go to the documentation of this file.
1 *> \brief \b DSYTRF_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 DSYTRF_AA_2STAGE + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrf_aa_2stage.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSYTRF_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 * DOUBLE PRECISION A( LDA, * ), TB( * ), WORK( * )
31 * ..
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> DSYTRF_AA_2STAGE computes the factorization of a real 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 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 DOUBLE PRECISION array, dimension (LDA,N)
70 *> On entry, the symmetric 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 DOUBLE PRECISION 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] WORK
106 *> \verbatim
107 *> WORK is DOUBLE PRECISION workspace of size LWORK
108 *> \endverbatim
109 *>
110 *> \param[in] LWORK
111 *> \verbatim
112 *> The size of WORK. LWORK >= N, internally used to select NB
113 *> such that LWORK >= N*NB.
114 *>
115 *> If LWORK = -1, then a workspace query is assumed; the
116 *> routine only calculates the optimal size of the WORK array,
117 *> returns this value as the first entry of the WORK array, and
118 *> no error message related to LWORK is issued by XERBLA.
119 *> \endverbatim
120 *>
121 *> \param[out] IPIV
122 *> \verbatim
123 *> IPIV is INTEGER array, dimension (N)
124 *> On exit, it contains the details of the interchanges, i.e.,
125 *> the row and column k of A were interchanged with the
126 *> row and column IPIV(k).
127 *> \endverbatim
128 *>
129 *> \param[out] IPIV2
130 *> \verbatim
131 *> IPIV is INTEGER array, dimension (N)
132 *> On exit, it contains the details of the interchanges, i.e.,
133 *> the row and column k of T were interchanged with the
134 *> row and column IPIV(k).
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 doubleSYcomputational
156 *
157 * =====================================================================
158  SUBROUTINE dsytrf_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  DOUBLE PRECISION A( lda, * ), TB( * ), WORK( * )
175 * ..
176 *
177 * =====================================================================
178 * .. Parameters ..
179  DOUBLE PRECISION ZERO, ONE
180  parameter( zero = 0.0d+0, one = 1.0d+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  DOUBLE PRECISION PIV
187 * ..
188 * .. External Functions ..
189  LOGICAL LSAME
190  INTEGER ILAENV
191  EXTERNAL lsame, ilaenv
192 * ..
193 * .. External Subroutines ..
194  EXTERNAL xerbla, dcopy, dlacgv, dlacpy,
195  $ dlaset, dgbtrf, dgemm, dgetrf,
196  $ dsygst, dswap, dtrsm
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( 'DSYTRF_AA_2STAGE', -info )
223  RETURN
224  END IF
225 *
226 * Answer the query
227 *
228  nb = ilaenv( 1, 'DSYTRF_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,I+1)*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 dgemm( '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 dgemm( '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 dlacpy( '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 dgemm( '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 dgemm( '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 dgemm( '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 dsygst( 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  DO k = i+1, kb
346  tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
347  $ = tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
348  END DO
349  END DO
350 *
351  IF( j.LT.nt-1 ) THEN
352  IF( j.GT.0 ) THEN
353 *
354 * Compute H(J,J)
355 *
356  IF( j.EQ.1 ) THEN
357  CALL dgemm( 'NoTranspose', 'NoTranspose',
358  $ kb, kb, kb,
359  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
360  $ a( (j-1)*nb+1, j*nb+1 ), lda,
361  $ zero, work( j*nb+1 ), n )
362  ELSE
363  CALL dgemm( 'NoTranspose', 'NoTranspose',
364  $ kb, kb, nb+kb,
365  $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
366  $ ldtb-1,
367  $ a( (j-2)*nb+1, j*nb+1 ), lda,
368  $ zero, work( j*nb+1 ), n )
369  END IF
370 *
371 * Update with the previous column
372 *
373  CALL dgemm( 'Transpose', 'NoTranspose',
374  $ nb, n-(j+1)*nb, j*nb,
375  $ -one, work( nb+1 ), n,
376  $ a( 1, (j+1)*nb+1 ), lda,
377  $ one, a( j*nb+1, (j+1)*nb+1 ), lda )
378  END IF
379 *
380 * Copy panel to workspace to call DGETRF
381 *
382  DO k = 1, nb
383  CALL dcopy( n-(j+1)*nb,
384  $ a( j*nb+k, (j+1)*nb+1 ), lda,
385  $ work( 1+(k-1)*n ), 1 )
386  END DO
387 *
388 * Factorize panel
389 *
390  CALL dgetrf( n-(j+1)*nb, nb,
391  $ work, n,
392  $ ipiv( (j+1)*nb+1 ), iinfo )
393 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
394 c INFO = IINFO+(J+1)*NB
395 c END IF
396 *
397 * Copy panel back
398 *
399  DO k = 1, nb
400  CALL dcopy( n-(j+1)*nb,
401  $ work( 1+(k-1)*n ), 1,
402  $ a( j*nb+k, (j+1)*nb+1 ), lda )
403  END DO
404 *
405 * Compute T(J+1, J), zero out for GEMM update
406 *
407  kb = min(nb, n-(j+1)*nb)
408  CALL dlaset( 'Full', kb, nb, zero, zero,
409  $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
410  CALL dlacpy( 'Upper', kb, nb,
411  $ work, n,
412  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
413  IF( j.GT.0 ) THEN
414  CALL dtrsm( 'R', 'U', 'N', 'U', kb, nb, one,
415  $ a( (j-1)*nb+1, j*nb+1 ), lda,
416  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
417  END IF
418 *
419 * Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
420 * updates
421 *
422  DO k = 1, nb
423  DO i = 1, kb
424  tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
425  $ = tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
426  END DO
427  END DO
428  CALL dlaset( 'Lower', kb, nb, zero, one,
429  $ a( j*nb+1, (j+1)*nb+1), lda )
430 *
431 * Apply pivots to trailing submatrix of A
432 *
433  DO k = 1, kb
434 * > Adjust ipiv
435  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
436 *
437  i1 = (j+1)*nb+k
438  i2 = ipiv( (j+1)*nb+k )
439  IF( i1.NE.i2 ) THEN
440 * > Apply pivots to previous columns of L
441  CALL dswap( k-1, a( (j+1)*nb+1, i1 ), 1,
442  $ a( (j+1)*nb+1, i2 ), 1 )
443 * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
444  CALL dswap( i2-i1-1, a( i1, i1+1 ), lda,
445  $ a( i1+1, i2 ), 1 )
446 * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
447  CALL dswap( n-i2, a( i1, i2+1 ), lda,
448  $ a( i2, i2+1 ), lda )
449 * > Swap A(I1, I1) with A(I2, I2)
450  piv = a( i1, i1 )
451  a( i1, i1 ) = a( i2, i2 )
452  a( i2, i2 ) = piv
453 * > Apply pivots to previous columns of L
454  IF( j.GT.0 ) THEN
455  CALL dswap( j*nb, a( 1, i1 ), 1,
456  $ a( 1, i2 ), 1 )
457  END IF
458  ENDIF
459  END DO
460  END IF
461  END DO
462  ELSE
463 *
464 * .....................................................
465 * Factorize A as L*D*L**T using the lower triangle of A
466 * .....................................................
467 *
468  DO j = 0, nt-1
469 *
470 * Generate Jth column of W and H
471 *
472  kb = min(nb, n-j*nb)
473  DO i = 1, j-1
474  IF( i.EQ.1 ) THEN
475 * H(I,J) = T(I,I)*L(J,I)' + T(I+1,I)'*L(J,I+1)'
476  IF( i .EQ. j-1) THEN
477  jb = nb+kb
478  ELSE
479  jb = 2*nb
480  END IF
481  CALL dgemm( 'NoTranspose', 'Transpose',
482  $ nb, kb, jb,
483  $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
484  $ a( j*nb+1, (i-1)*nb+1 ), lda,
485  $ zero, work( i*nb+1 ), n )
486  ELSE
487 * 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)'
488  IF( i .EQ. j-1) THEN
489  jb = 2*nb+kb
490  ELSE
491  jb = 3*nb
492  END IF
493  CALL dgemm( 'NoTranspose', 'Transpose',
494  $ nb, kb, jb,
495  $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
496  $ ldtb-1,
497  $ a( j*nb+1, (i-2)*nb+1 ), lda,
498  $ zero, work( i*nb+1 ), n )
499  END IF
500  END DO
501 *
502 * Compute T(J,J)
503 *
504  CALL dlacpy( 'Lower', kb, kb, a( j*nb+1, j*nb+1 ), lda,
505  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
506  IF( j.GT.1 ) THEN
507 * T(J,J) = L(J,1:J)*H(1:J)
508  CALL dgemm( 'NoTranspose', 'NoTranspose',
509  $ kb, kb, (j-1)*nb,
510  $ -one, a( j*nb+1, 1 ), lda,
511  $ work( nb+1 ), n,
512  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
513 * T(J,J) += L(J,J)*T(J,J-1)*L(J,J-1)'
514  CALL dgemm( 'NoTranspose', 'NoTranspose',
515  $ kb, nb, kb,
516  $ one, a( j*nb+1, (j-1)*nb+1 ), lda,
517  $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
518  $ zero, work( 1 ), n )
519  CALL dgemm( 'NoTranspose', 'Transpose',
520  $ kb, kb, nb,
521  $ -one, work( 1 ), n,
522  $ a( j*nb+1, (j-2)*nb+1 ), lda,
523  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
524  END IF
525  IF( j.GT.0 ) THEN
526  CALL dsygst( 1, 'Lower', kb,
527  $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
528  $ a( j*nb+1, (j-1)*nb+1 ), lda, iinfo )
529  END IF
530 *
531 * Expand T(J,J) into full format
532 *
533  DO i = 1, kb
534  DO k = i+1, kb
535  tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
536  $ = tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
537  END DO
538  END DO
539 *
540  IF( j.LT.nt-1 ) THEN
541  IF( j.GT.0 ) THEN
542 *
543 * Compute H(J,J)
544 *
545  IF( j.EQ.1 ) THEN
546  CALL dgemm( 'NoTranspose', 'Transpose',
547  $ kb, kb, kb,
548  $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
549  $ a( j*nb+1, (j-1)*nb+1 ), lda,
550  $ zero, work( j*nb+1 ), n )
551  ELSE
552  CALL dgemm( 'NoTranspose', 'Transpose',
553  $ kb, kb, nb+kb,
554  $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
555  $ ldtb-1,
556  $ a( j*nb+1, (j-2)*nb+1 ), lda,
557  $ zero, work( j*nb+1 ), n )
558  END IF
559 *
560 * Update with the previous column
561 *
562  CALL dgemm( 'NoTranspose', 'NoTranspose',
563  $ n-(j+1)*nb, nb, j*nb,
564  $ -one, a( (j+1)*nb+1, 1 ), lda,
565  $ work( nb+1 ), n,
566  $ one, a( (j+1)*nb+1, j*nb+1 ), lda )
567  END IF
568 *
569 * Factorize panel
570 *
571  CALL dgetrf( n-(j+1)*nb, nb,
572  $ a( (j+1)*nb+1, j*nb+1 ), lda,
573  $ ipiv( (j+1)*nb+1 ), iinfo )
574 c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
575 c INFO = IINFO+(J+1)*NB
576 c END IF
577 *
578 * Compute T(J+1, J), zero out for GEMM update
579 *
580  kb = min(nb, n-(j+1)*nb)
581  CALL dlaset( 'Full', kb, nb, zero, zero,
582  $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
583  CALL dlacpy( 'Upper', kb, nb,
584  $ a( (j+1)*nb+1, j*nb+1 ), lda,
585  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
586  IF( j.GT.0 ) THEN
587  CALL dtrsm( 'R', 'L', 'T', 'U', kb, nb, one,
588  $ a( j*nb+1, (j-1)*nb+1 ), lda,
589  $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
590  END IF
591 *
592 * Copy T(J+1,J) into T(J, J+1), both upper/lower for GEMM
593 * updates
594 *
595  DO k = 1, nb
596  DO i = 1, kb
597  tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
598  $ = tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
599  END DO
600  END DO
601  CALL dlaset( 'Upper', kb, nb, zero, one,
602  $ a( (j+1)*nb+1, j*nb+1), lda )
603 *
604 * Apply pivots to trailing submatrix of A
605 *
606  DO k = 1, kb
607 * > Adjust ipiv
608  ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
609 *
610  i1 = (j+1)*nb+k
611  i2 = ipiv( (j+1)*nb+k )
612  IF( i1.NE.i2 ) THEN
613 * > Apply pivots to previous columns of L
614  CALL dswap( k-1, a( i1, (j+1)*nb+1 ), lda,
615  $ a( i2, (j+1)*nb+1 ), lda )
616 * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
617  CALL dswap( i2-i1-1, a( i1+1, i1 ), 1,
618  $ a( i2, i1+1 ), lda )
619 * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
620  CALL dswap( n-i2, a( i2+1, i1 ), 1,
621  $ a( i2+1, i2 ), 1 )
622 * > Swap A(I1, I1) with A(I2, I2)
623  piv = a( i1, i1 )
624  a( i1, i1 ) = a( i2, i2 )
625  a( i2, i2 ) = piv
626 * > Apply pivots to previous columns of L
627  IF( j.GT.0 ) THEN
628  CALL dswap( j*nb, a( i1, 1 ), lda,
629  $ a( i2, 1 ), lda )
630  END IF
631  ENDIF
632  END DO
633 *
634 * Apply pivots to previous columns of L
635 *
636 c CALL DLASWP( J*NB, A( 1, 1 ), LDA,
637 c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
638  END IF
639  END DO
640  END IF
641 *
642 * Factor the band matrix
643  CALL dgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
644 *
645 * End of DSYTRF_AA_2STAGE
646 *
647  END
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dgetrf(M, N, A, LDA, IPIV, INFO)
DGETRF
Definition: dgetrf.f:110
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:183
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
subroutine dsytrf_aa_2stage(UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2, WORK, LWORK, INFO)
DSYTRF_AA_2STAGE
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
subroutine dgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
DGBTRF
Definition: dgbtrf.f:146
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dsygst(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
DSYGST
Definition: dsygst.f:129