LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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,
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 )
392c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
393c INFO = IINFO+(J+1)*NB
394c 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 )
575c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
576c INFO = IINFO+(J+1)*NB
577c 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*
639c CALL SLASWP( J*NB, A( 1, 1 ), LDA,
640c $ (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