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