LAPACK 3.12.1
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*> Download SSYTRF_AA_2STAGE + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytrf_aa_2stage.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytrf_aa_2stage.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytrf_aa_2stage.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
20* IPIV2, WORK, LWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER UPLO
24* INTEGER N, LDA, LTB, LWORK, INFO
25* ..
26* .. Array Arguments ..
27* INTEGER IPIV( * ), IPIV2( * )
28* REAL A( LDA, * ), TB( * ), WORK( * )
29* ..
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> SSYTRF_AA_2STAGE computes the factorization of a real symmetric matrix A
37*> using the Aasen's algorithm. The form of the factorization is
38*>
39*> A = U**T*T*U or A = L*T*L**T
40*>
41*> where U (or L) is a product of permutation and unit upper (lower)
42*> triangular matrices, and T is a symmetric band matrix with the
43*> bandwidth of NB (NB is internally selected and stored in TB( 1 ), and T is
44*> LU factorized with partial pivoting).
45*>
46*> This is the blocked version of the algorithm, calling Level 3 BLAS.
47*> \endverbatim
48*
49* Arguments:
50* ==========
51*
52*> \param[in] UPLO
53*> \verbatim
54*> UPLO is CHARACTER*1
55*> = 'U': Upper triangle of A is stored;
56*> = 'L': Lower triangle of A is stored.
57*> \endverbatim
58*>
59*> \param[in] N
60*> \verbatim
61*> N is INTEGER
62*> The order of the matrix A. N >= 0.
63*> \endverbatim
64*>
65*> \param[in,out] A
66*> \verbatim
67*> A is REAL array, dimension (LDA,N)
68*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
69*> N-by-N upper triangular part of A contains the upper
70*> triangular part of the matrix A, and the strictly lower
71*> triangular part of A is not referenced. If UPLO = 'L', the
72*> leading N-by-N lower triangular part of A contains the lower
73*> triangular part of the matrix A, and the strictly upper
74*> triangular part of A is not referenced.
75*>
76*> On exit, L is stored below (or above) the subdiagonal blocks,
77*> when UPLO is 'L' (or 'U').
78*> \endverbatim
79*>
80*> \param[in] LDA
81*> \verbatim
82*> LDA is INTEGER
83*> The leading dimension of the array A. LDA >= max(1,N).
84*> \endverbatim
85*>
86*> \param[out] TB
87*> \verbatim
88*> TB is REAL array, dimension (LTB)
89*> On exit, details of the LU factorization of the band matrix.
90*> \endverbatim
91*>
92*> \param[in] LTB
93*> \verbatim
94*> LTB is INTEGER
95*> The size of the array TB. LTB >= MAX(1,4*N), internally
96*> used to select NB such that LTB >= (3*NB+1)*N.
97*>
98*> If LTB = -1, then a workspace query is assumed; the
99*> routine only calculates the optimal size of LTB,
100*> returns this value as the first entry of TB, and
101*> no error message related to LTB is issued by XERBLA.
102*> \endverbatim
103*>
104*> \param[out] IPIV
105*> \verbatim
106*> IPIV is INTEGER array, dimension (N)
107*> On exit, it contains the details of the interchanges, i.e.,
108*> the row and column k of A were interchanged with the
109*> row and column IPIV(k).
110*> \endverbatim
111*>
112*> \param[out] IPIV2
113*> \verbatim
114*> IPIV2 is INTEGER array, dimension (N)
115*> On exit, it contains the details of the interchanges, i.e.,
116*> the row and column k of T were interchanged with the
117*> row and column IPIV(k).
118*> \endverbatim
119*>
120*> \param[out] WORK
121*> \verbatim
122*> WORK is REAL workspace of size (MAX(1,LWORK))
123*> \endverbatim
124*>
125*> \param[in] LWORK
126*> \verbatim
127*> LWORK is INTEGER
128*> The size of WORK. LWORK >= MAX(1,N), internally used to
129*> select NB 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*> \ingroup hetrf_aa_2stage
154*
155* =====================================================================
156 SUBROUTINE ssytrf_aa_2stage( UPLO, N, A, LDA, TB, LTB, IPIV,
157 $ IPIV2, WORK, LWORK, INFO )
158*
159* -- LAPACK computational routine --
160* -- LAPACK is a software package provided by Univ. of Tennessee, --
161* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162*
163 IMPLICIT NONE
164*
165* .. Scalar Arguments ..
166 CHARACTER UPLO
167 INTEGER N, LDA, LTB, LWORK, INFO
168* ..
169* .. Array Arguments ..
170 INTEGER IPIV( * ), IPIV2( * )
171 REAL A( LDA, * ), TB( * ), WORK( * )
172* ..
173*
174* =====================================================================
175* .. Parameters ..
176 REAL ZERO, ONE
177 parameter( zero = 0.0e+0, one = 1.0e+0 )
178*
179* .. Local Scalars ..
180 LOGICAL UPPER, TQUERY, WQUERY
181 INTEGER I, J, K, I1, I2, TD
182 INTEGER LDTB, NB, KB, JB, NT, IINFO
183 REAL PIV
184* ..
185* .. External Functions ..
186 LOGICAL LSAME
187 INTEGER ILAENV
188 REAL SROUNDUP_LWORK
189 EXTERNAL lsame, ilaenv, sroundup_lwork
190* ..
191* .. External Subroutines ..
192 EXTERNAL xerbla, scopy, slacpy,
194 $ ssygst, sswap, strsm
195* ..
196* .. Intrinsic Functions ..
197 INTRINSIC min, max
198* ..
199* .. Executable Statements ..
200*
201* Test the input parameters.
202*
203 info = 0
204 upper = lsame( uplo, 'U' )
205 wquery = ( lwork.EQ.-1 )
206 tquery = ( ltb.EQ.-1 )
207 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
208 info = -1
209 ELSE IF( n.LT.0 ) THEN
210 info = -2
211 ELSE IF( lda.LT.max( 1, n ) ) THEN
212 info = -4
213 ELSE IF( ltb.LT.max( 1, 4*n ) .AND. .NOT.tquery ) THEN
214 info = -6
215 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.wquery ) THEN
216 info = -10
217 END IF
218*
219 IF( info.NE.0 ) THEN
220 CALL xerbla( 'SSYTRF_AA_2STAGE', -info )
221 RETURN
222 END IF
223*
224* Answer the query
225*
226 nb = ilaenv( 1, 'SSYTRF_AA_2STAGE', uplo, n, -1, -1, -1 )
227 IF( info.EQ.0 ) THEN
228 IF( tquery ) THEN
229 tb( 1 ) = sroundup_lwork( max( 1, (3*nb+1)*n ) )
230 END IF
231 IF( wquery ) THEN
232 work( 1 ) = sroundup_lwork( max( 1, n*nb ) )
233 END IF
234 END IF
235 IF( tquery .OR. wquery ) THEN
236 RETURN
237 END IF
238*
239* Quick return
240*
241 IF( n.EQ.0 ) THEN
242 RETURN
243 ENDIF
244*
245* Determine the number of the block size
246*
247 ldtb = ltb/n
248 IF( ldtb .LT. 3*nb+1 ) THEN
249 nb = (ldtb-1)/3
250 END IF
251 IF( lwork .LT. nb*n ) THEN
252 nb = lwork/n
253 END IF
254*
255* Determine the number of the block columns
256*
257 nt = (n+nb-1)/nb
258 td = 2*nb
259 kb = min(nb, n)
260*
261* Initialize vectors/matrices
262*
263 DO j = 1, kb
264 ipiv( j ) = j
265 END DO
266*
267* Save NB
268*
269 tb( 1 ) = real( nb )
270*
271 IF( upper ) THEN
272*
273* .....................................................
274* Factorize A as U**T*D*U using the upper triangle of A
275* .....................................................
276*
277 DO j = 0, nt-1
278*
279* Generate Jth column of W and H
280*
281 kb = min(nb, n-j*nb)
282 DO i = 1, j-1
283 IF( i.EQ.1 ) THEN
284* H(I,J) = T(I,I)*U(I,J) + T(I+1,I)*U(I+1,J)
285 IF( i .EQ. (j-1) ) THEN
286 jb = nb+kb
287 ELSE
288 jb = 2*nb
289 END IF
290 CALL sgemm( 'NoTranspose', 'NoTranspose',
291 $ nb, kb, jb,
292 $ one, tb( td+1 + (i*nb)*ldtb ), ldtb-1,
293 $ a( (i-1)*nb+1, j*nb+1 ), lda,
294 $ zero, work( i*nb+1 ), n )
295 ELSE
296* 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)
297 IF( i .EQ. j-1) THEN
298 jb = 2*nb+kb
299 ELSE
300 jb = 3*nb
301 END IF
302 CALL sgemm( 'NoTranspose', 'NoTranspose',
303 $ nb, kb, jb,
304 $ one, tb( td+nb+1 + ((i-1)*nb)*ldtb ),
305 $ ldtb-1,
306 $ a( (i-2)*nb+1, j*nb+1 ), lda,
307 $ zero, work( i*nb+1 ), n )
308 END IF
309 END DO
310*
311* Compute T(J,J)
312*
313 CALL slacpy( 'Upper', kb, kb, a( j*nb+1, j*nb+1 ), lda,
314 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
315 IF( j.GT.1 ) THEN
316* T(J,J) = U(1:J,J)'*H(1:J)
317 CALL sgemm( 'Transpose', 'NoTranspose',
318 $ kb, kb, (j-1)*nb,
319 $ -one, a( 1, j*nb+1 ), lda,
320 $ work( nb+1 ), n,
321 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
322* T(J,J) += U(J,J)'*T(J,J-1)*U(J-1,J)
323 CALL sgemm( 'Transpose', 'NoTranspose',
324 $ kb, nb, kb,
325 $ one, a( (j-1)*nb+1, j*nb+1 ), lda,
326 $ tb( td+nb+1 + ((j-1)*nb)*ldtb ), ldtb-1,
327 $ zero, work( 1 ), n )
328 CALL sgemm( 'NoTranspose', 'NoTranspose',
329 $ kb, kb, nb,
330 $ -one, work( 1 ), n,
331 $ a( (j-2)*nb+1, j*nb+1 ), lda,
332 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1 )
333 END IF
334 IF( j.GT.0 ) THEN
335 CALL ssygst( 1, 'Upper', kb,
336 $ tb( td+1 + (j*nb)*ldtb ), ldtb-1,
337 $ a( (j-1)*nb+1, j*nb+1 ), lda, iinfo )
338 END IF
339*
340* Expand T(J,J) into full format
341*
342 DO i = 1, kb
343 DO k = i+1, kb
344 tb( td+(k-i)+1 + (j*nb+i-1)*ldtb )
345 $ = tb( td-(k-(i+1)) + (j*nb+k-1)*ldtb )
346 END DO
347 END DO
348*
349 IF( j.LT.nt-1 ) THEN
350 IF( j.GT.0 ) THEN
351*
352* Compute H(J,J)
353*
354 IF( j.EQ.1 ) THEN
355 CALL sgemm( 'NoTranspose', 'NoTranspose',
356 $ kb, kb, kb,
357 $ one, tb( td+1 + (j*nb)*ldtb ), ldtb-1,
358 $ a( (j-1)*nb+1, j*nb+1 ), lda,
359 $ zero, work( j*nb+1 ), n )
360 ELSE
361 CALL sgemm( 'NoTranspose', 'NoTranspose',
362 $ kb, kb, nb+kb,
363 $ one, tb( td+nb+1 + ((j-1)*nb)*ldtb ),
364 $ ldtb-1,
365 $ a( (j-2)*nb+1, j*nb+1 ), lda,
366 $ zero, work( j*nb+1 ), n )
367 END IF
368*
369* Update with the previous column
370*
371 CALL sgemm( 'Transpose', 'NoTranspose',
372 $ nb, n-(j+1)*nb, j*nb,
373 $ -one, work( nb+1 ), n,
374 $ a( 1, (j+1)*nb+1 ), lda,
375 $ one, a( j*nb+1, (j+1)*nb+1 ), lda )
376 END IF
377*
378* Copy panel to workspace to call SGETRF
379*
380 DO k = 1, nb
381 CALL scopy( n-(j+1)*nb,
382 $ a( j*nb+k, (j+1)*nb+1 ), lda,
383 $ work( 1+(k-1)*n ), 1 )
384 END DO
385*
386* Factorize panel
387*
388 CALL sgetrf( n-(j+1)*nb, nb,
389 $ work, n,
390 $ ipiv( (j+1)*nb+1 ), iinfo )
391c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
392c INFO = IINFO+(J+1)*NB
393c END IF
394*
395* Copy panel back
396*
397 DO k = 1, nb
398 CALL scopy( n-(j+1)*nb,
399 $ work( 1+(k-1)*n ), 1,
400 $ a( j*nb+k, (j+1)*nb+1 ), lda )
401 END DO
402*
403* Compute T(J+1, J), zero out for GEMM update
404*
405 kb = min(nb, n-(j+1)*nb)
406 CALL slaset( 'Full', kb, nb, zero, zero,
407 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
408 CALL slacpy( 'Upper', kb, nb,
409 $ work, n,
410 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
411 IF( j.GT.0 ) THEN
412 CALL strsm( 'R', 'U', 'N', 'U', kb, nb, one,
413 $ a( (j-1)*nb+1, j*nb+1 ), lda,
414 $ tb( td+nb+1 + (j*nb)*ldtb ), ldtb-1 )
415 END IF
416*
417* Copy T(J,J+1) into T(J+1, J), both upper/lower for GEMM
418* updates
419*
420 DO k = 1, nb
421 DO i = 1, kb
422 tb( td-nb+k-i+1 + (j*nb+nb+i-1)*ldtb )
423 $ = tb( td+nb+i-k+1 + (j*nb+k-1)*ldtb )
424 END DO
425 END DO
426 CALL slaset( 'Lower', kb, nb, zero, one,
427 $ a( j*nb+1, (j+1)*nb+1), lda )
428*
429* Apply pivots to trailing submatrix of A
430*
431 DO k = 1, kb
432* > Adjust ipiv
433 ipiv( (j+1)*nb+k ) = ipiv( (j+1)*nb+k ) + (j+1)*nb
434*
435 i1 = (j+1)*nb+k
436 i2 = ipiv( (j+1)*nb+k )
437 IF( i1.NE.i2 ) THEN
438* > Apply pivots to previous columns of L
439 CALL sswap( k-1, a( (j+1)*nb+1, i1 ), 1,
440 $ a( (j+1)*nb+1, i2 ), 1 )
441* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
442 IF( i2.GT.(i1+1) )
443 $ CALL sswap( i2-i1-1, a( i1, i1+1 ), lda,
444 $ a( i1+1, i2 ), 1 )
445* > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
446 IF( i2.LT.n )
447 $ CALL sswap( 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 sswap( 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 sgemm( '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 sgemm( '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 slacpy( '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 sgemm( '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 sgemm( '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 sgemm( '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 ssygst( 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 sgemm( '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 sgemm( '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 sgemm( '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 sgetrf( n-(j+1)*nb, nb,
572 $ a( (j+1)*nb+1, j*nb+1 ), lda,
573 $ ipiv( (j+1)*nb+1 ), iinfo )
574c IF (IINFO.NE.0 .AND. INFO.EQ.0) THEN
575c INFO = IINFO+(J+1)*NB
576c 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 slaset( 'Full', kb, nb, zero, zero,
582 $ tb( td+nb+1 + (j*nb)*ldtb), ldtb-1 )
583 CALL slacpy( '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 strsm( '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 slaset( '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 sswap( 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 IF( i2.GT.(i1+1) )
618 $ CALL sswap( i2-i1-1, a( i1+1, i1 ), 1,
619 $ a( i2, i1+1 ), lda )
620* > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
621 IF( i2.LT.n )
622 $ CALL sswap( n-i2, a( i2+1, i1 ), 1,
623 $ a( i2+1, i2 ), 1 )
624* > Swap A(I1, I1) with A(I2, I2)
625 piv = a( i1, i1 )
626 a( i1, i1 ) = a( i2, i2 )
627 a( i2, i2 ) = piv
628* > Apply pivots to previous columns of L
629 IF( j.GT.0 ) THEN
630 CALL sswap( j*nb, a( i1, 1 ), lda,
631 $ a( i2, 1 ), lda )
632 END IF
633 ENDIF
634 END DO
635*
636* Apply pivots to previous columns of L
637*
638c CALL SLASWP( J*NB, A( 1, 1 ), LDA,
639c $ (J+1)*NB+1, (J+1)*NB+KB, IPIV, 1 )
640 END IF
641 END DO
642 END IF
643*
644* Factor the band matrix
645 CALL sgbtrf( n, n, nb, nb, tb, ldtb, ipiv2, info )
646*
647 RETURN
648*
649* End of SSYTRF_AA_2STAGE
650*
651 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgbtrf(m, n, kl, ku, ab, ldab, ipiv, info)
SGBTRF
Definition sgbtrf.f:142
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine sgetrf(m, n, a, lda, ipiv, info)
SGETRF
Definition sgetrf.f:106
subroutine ssygst(itype, uplo, n, a, lda, b, ldb, info)
SSYGST
Definition ssygst.f:125
subroutine ssytrf_aa_2stage(uplo, n, a, lda, tb, ltb, ipiv, ipiv2, work, lwork, info)
SSYTRF_AA_2STAGE
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
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:108
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181