LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slasyf_rook.f
Go to the documentation of this file.
1*> \brief \b SLASYF_ROOK computes a partial factorization of a real symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLASYF_ROOK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasyf_rook.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasyf_rook.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasyf_rook.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, KB, LDA, LDW, N, NB
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* REAL A( LDA, * ), W( LDW, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SLASYF_ROOK computes a partial factorization of a real symmetric
39*> matrix A using the bounded Bunch-Kaufman ("rook") diagonal
40*> pivoting method. The partial factorization has the form:
41*>
42*> A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
43*> ( 0 U22 ) ( 0 D ) ( U12**T U22**T )
44*>
45*> A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) if UPLO = 'L'
46*> ( L21 I ) ( 0 A22 ) ( 0 I )
47*>
48*> where the order of D is at most NB. The actual order is returned in
49*> the argument KB, and is either NB or NB-1, or N if N <= NB.
50*>
51*> SLASYF_ROOK is an auxiliary routine called by SSYTRF_ROOK. It uses
52*> blocked code (calling Level 3 BLAS) to update the submatrix
53*> A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
54*> \endverbatim
55*
56* Arguments:
57* ==========
58*
59*> \param[in] UPLO
60*> \verbatim
61*> UPLO is CHARACTER*1
62*> Specifies whether the upper or lower triangular part of the
63*> symmetric matrix A is stored:
64*> = 'U': Upper triangular
65*> = 'L': Lower triangular
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*> N is INTEGER
71*> The order of the matrix A. N >= 0.
72*> \endverbatim
73*>
74*> \param[in] NB
75*> \verbatim
76*> NB is INTEGER
77*> The maximum number of columns of the matrix A that should be
78*> factored. NB should be at least 2 to allow for 2-by-2 pivot
79*> blocks.
80*> \endverbatim
81*>
82*> \param[out] KB
83*> \verbatim
84*> KB is INTEGER
85*> The number of columns of A that were actually factored.
86*> KB is either NB-1 or NB, or N if N <= NB.
87*> \endverbatim
88*>
89*> \param[in,out] A
90*> \verbatim
91*> A is REAL array, dimension (LDA,N)
92*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
93*> n-by-n upper triangular part of A contains the upper
94*> triangular part of the matrix A, and the strictly lower
95*> triangular part of A is not referenced. If UPLO = 'L', the
96*> leading n-by-n lower triangular part of A contains the lower
97*> triangular part of the matrix A, and the strictly upper
98*> triangular part of A is not referenced.
99*> On exit, A contains details of the partial factorization.
100*> \endverbatim
101*>
102*> \param[in] LDA
103*> \verbatim
104*> LDA is INTEGER
105*> The leading dimension of the array A. LDA >= max(1,N).
106*> \endverbatim
107*>
108*> \param[out] IPIV
109*> \verbatim
110*> IPIV is INTEGER array, dimension (N)
111*> Details of the interchanges and the block structure of D.
112*>
113*> If UPLO = 'U':
114*> Only the last KB elements of IPIV are set.
115*>
116*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
117*> interchanged and D(k,k) is a 1-by-1 diagonal block.
118*>
119*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
120*> columns k and -IPIV(k) were interchanged and rows and
121*> columns k-1 and -IPIV(k-1) were inerchaged,
122*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
123*>
124*> If UPLO = 'L':
125*> Only the first KB elements of IPIV are set.
126*>
127*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
128*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
129*>
130*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
131*> columns k and -IPIV(k) were interchanged and rows and
132*> columns k+1 and -IPIV(k+1) were inerchaged,
133*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
134*> \endverbatim
135*>
136*> \param[out] W
137*> \verbatim
138*> W is REAL array, dimension (LDW,NB)
139*> \endverbatim
140*>
141*> \param[in] LDW
142*> \verbatim
143*> LDW is INTEGER
144*> The leading dimension of the array W. LDW >= max(1,N).
145*> \endverbatim
146*>
147*> \param[out] INFO
148*> \verbatim
149*> INFO is INTEGER
150*> = 0: successful exit
151*> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
152*> has been completed, but the block diagonal matrix D is
153*> exactly singular.
154*> \endverbatim
155*
156* Authors:
157* ========
158*
159*> \author Univ. of Tennessee
160*> \author Univ. of California Berkeley
161*> \author Univ. of Colorado Denver
162*> \author NAG Ltd.
163*
164*> \ingroup realSYcomputational
165*
166*> \par Contributors:
167* ==================
168*>
169*> \verbatim
170*>
171*> November 2013, Igor Kozachenko,
172*> Computer Science Division,
173*> University of California, Berkeley
174*>
175*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
176*> School of Mathematics,
177*> University of Manchester
178*>
179*> \endverbatim
180*
181* =====================================================================
182 SUBROUTINE slasyf_rook( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
183 $ INFO )
184*
185* -- LAPACK computational routine --
186* -- LAPACK is a software package provided by Univ. of Tennessee, --
187* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188*
189* .. Scalar Arguments ..
190 CHARACTER UPLO
191 INTEGER INFO, KB, LDA, LDW, N, NB
192* ..
193* .. Array Arguments ..
194 INTEGER IPIV( * )
195 REAL A( LDA, * ), W( LDW, * )
196* ..
197*
198* =====================================================================
199*
200* .. Parameters ..
201 REAL ZERO, ONE
202 parameter( zero = 0.0e+0, one = 1.0e+0 )
203 REAL EIGHT, SEVTEN
204 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
205* ..
206* .. Local Scalars ..
207 LOGICAL DONE
208 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK,
209 $ kw, kkw, kp, kstep, p, ii
210
211 REAL ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
212 $ stemp, r1, rowmax, t, sfmin
213* ..
214* .. External Functions ..
215 LOGICAL LSAME
216 INTEGER ISAMAX
217 REAL SLAMCH
218 EXTERNAL lsame, isamax, slamch
219* ..
220* .. External Subroutines ..
221 EXTERNAL scopy, sgemm, sgemv, sscal, sswap
222* ..
223* .. Intrinsic Functions ..
224 INTRINSIC abs, max, min, sqrt
225* ..
226* .. Executable Statements ..
227*
228 info = 0
229*
230* Initialize ALPHA for use in choosing pivot block size.
231*
232 alpha = ( one+sqrt( sevten ) ) / eight
233*
234* Compute machine safe minimum
235*
236 sfmin = slamch( 'S' )
237*
238 IF( lsame( uplo, 'U' ) ) THEN
239*
240* Factorize the trailing columns of A using the upper triangle
241* of A and working backwards, and compute the matrix W = U12*D
242* for use in updating A11
243*
244* K is the main loop index, decreasing from N in steps of 1 or 2
245*
246 k = n
247 10 CONTINUE
248*
249* KW is the column of W which corresponds to column K of A
250*
251 kw = nb + k - n
252*
253* Exit from loop
254*
255 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
256 $ GO TO 30
257*
258 kstep = 1
259 p = k
260*
261* Copy column K of A to column KW of W and update it
262*
263 CALL scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
264 IF( k.LT.n )
265 $ CALL sgemv( 'No transpose', k, n-k, -one, a( 1, k+1 ),
266 $ lda, w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
267*
268* Determine rows and columns to be interchanged and whether
269* a 1-by-1 or 2-by-2 pivot block will be used
270*
271 absakk = abs( w( k, kw ) )
272*
273* IMAX is the row-index of the largest off-diagonal element in
274* column K, and COLMAX is its absolute value.
275* Determine both COLMAX and IMAX.
276*
277 IF( k.GT.1 ) THEN
278 imax = isamax( k-1, w( 1, kw ), 1 )
279 colmax = abs( w( imax, kw ) )
280 ELSE
281 colmax = zero
282 END IF
283*
284 IF( max( absakk, colmax ).EQ.zero ) THEN
285*
286* Column K is zero or underflow: set INFO and continue
287*
288 IF( info.EQ.0 )
289 $ info = k
290 kp = k
291 CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
292 ELSE
293*
294* ============================================================
295*
296* Test for interchange
297*
298* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
299* (used to handle NaN and Inf)
300*
301 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
302*
303* no interchange, use 1-by-1 pivot block
304*
305 kp = k
306*
307 ELSE
308*
309 done = .false.
310*
311* Loop until pivot found
312*
313 12 CONTINUE
314*
315* Begin pivot search loop body
316*
317*
318* Copy column IMAX to column KW-1 of W and update it
319*
320 CALL scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
321 CALL scopy( k-imax, a( imax, imax+1 ), lda,
322 $ w( imax+1, kw-1 ), 1 )
323*
324 IF( k.LT.n )
325 $ CALL sgemv( 'No transpose', k, n-k, -one,
326 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
327 $ one, w( 1, kw-1 ), 1 )
328*
329* JMAX is the column-index of the largest off-diagonal
330* element in row IMAX, and ROWMAX is its absolute value.
331* Determine both ROWMAX and JMAX.
332*
333 IF( imax.NE.k ) THEN
334 jmax = imax + isamax( k-imax, w( imax+1, kw-1 ),
335 $ 1 )
336 rowmax = abs( w( jmax, kw-1 ) )
337 ELSE
338 rowmax = zero
339 END IF
340*
341 IF( imax.GT.1 ) THEN
342 itemp = isamax( imax-1, w( 1, kw-1 ), 1 )
343 stemp = abs( w( itemp, kw-1 ) )
344 IF( stemp.GT.rowmax ) THEN
345 rowmax = stemp
346 jmax = itemp
347 END IF
348 END IF
349*
350* Equivalent to testing for
351* ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
352* (used to handle NaN and Inf)
353*
354 IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
355 $ THEN
356*
357* interchange rows and columns K and IMAX,
358* use 1-by-1 pivot block
359*
360 kp = imax
361*
362* copy column KW-1 of W to column KW of W
363*
364 CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
365*
366 done = .true.
367*
368* Equivalent to testing for ROWMAX.EQ.COLMAX,
369* (used to handle NaN and Inf)
370*
371 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
372 $ THEN
373*
374* interchange rows and columns K-1 and IMAX,
375* use 2-by-2 pivot block
376*
377 kp = imax
378 kstep = 2
379 done = .true.
380 ELSE
381*
382* Pivot not found: set params and repeat
383*
384 p = imax
385 colmax = rowmax
386 imax = jmax
387*
388* Copy updated JMAXth (next IMAXth) column to Kth of W
389*
390 CALL scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
391*
392 END IF
393*
394* End pivot search loop body
395*
396 IF( .NOT. done ) GOTO 12
397*
398 END IF
399*
400* ============================================================
401*
402 kk = k - kstep + 1
403*
404* KKW is the column of W which corresponds to column KK of A
405*
406 kkw = nb + kk - n
407*
408 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
409*
410* Copy non-updated column K to column P
411*
412 CALL scopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
413 CALL scopy( p, a( 1, k ), 1, a( 1, p ), 1 )
414*
415* Interchange rows K and P in last N-K+1 columns of A
416* and last N-K+2 columns of W
417*
418 CALL sswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
419 CALL sswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
420 END IF
421*
422* Updated column KP is already stored in column KKW of W
423*
424 IF( kp.NE.kk ) THEN
425*
426* Copy non-updated column KK to column KP
427*
428 a( kp, k ) = a( kk, k )
429 CALL scopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
430 $ lda )
431 CALL scopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
432*
433* Interchange rows KK and KP in last N-KK+1 columns
434* of A and W
435*
436 CALL sswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
437 CALL sswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
438 $ ldw )
439 END IF
440*
441 IF( kstep.EQ.1 ) THEN
442*
443* 1-by-1 pivot block D(k): column KW of W now holds
444*
445* W(k) = U(k)*D(k)
446*
447* where U(k) is the k-th column of U
448*
449* Store U(k) in column k of A
450*
451 CALL scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
452 IF( k.GT.1 ) THEN
453 IF( abs( a( k, k ) ).GE.sfmin ) THEN
454 r1 = one / a( k, k )
455 CALL sscal( k-1, r1, a( 1, k ), 1 )
456 ELSE IF( a( k, k ).NE.zero ) THEN
457 DO 14 ii = 1, k - 1
458 a( ii, k ) = a( ii, k ) / a( k, k )
459 14 CONTINUE
460 END IF
461 END IF
462*
463 ELSE
464*
465* 2-by-2 pivot block D(k): columns KW and KW-1 of W now
466* hold
467*
468* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
469*
470* where U(k) and U(k-1) are the k-th and (k-1)-th columns
471* of U
472*
473 IF( k.GT.2 ) THEN
474*
475* Store U(k) and U(k-1) in columns k and k-1 of A
476*
477 d12 = w( k-1, kw )
478 d11 = w( k, kw ) / d12
479 d22 = w( k-1, kw-1 ) / d12
480 t = one / ( d11*d22-one )
481 DO 20 j = 1, k - 2
482 a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
483 $ d12 )
484 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
485 $ d12 )
486 20 CONTINUE
487 END IF
488*
489* Copy D(k) to A
490*
491 a( k-1, k-1 ) = w( k-1, kw-1 )
492 a( k-1, k ) = w( k-1, kw )
493 a( k, k ) = w( k, kw )
494 END IF
495 END IF
496*
497* Store details of the interchanges in IPIV
498*
499 IF( kstep.EQ.1 ) THEN
500 ipiv( k ) = kp
501 ELSE
502 ipiv( k ) = -p
503 ipiv( k-1 ) = -kp
504 END IF
505*
506* Decrease K and return to the start of the main loop
507*
508 k = k - kstep
509 GO TO 10
510*
511 30 CONTINUE
512*
513* Update the upper triangle of A11 (= A(1:k,1:k)) as
514*
515* A11 := A11 - U12*D*U12**T = A11 - U12*W**T
516*
517* computing blocks of NB columns at a time
518*
519 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
520 jb = min( nb, k-j+1 )
521*
522* Update the upper triangle of the diagonal block
523*
524 DO 40 jj = j, j + jb - 1
525 CALL sgemv( 'No transpose', jj-j+1, n-k, -one,
526 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
527 $ a( j, jj ), 1 )
528 40 CONTINUE
529*
530* Update the rectangular superdiagonal block
531*
532 IF( j.GE.2 )
533 $ CALL sgemm( 'No transpose', 'Transpose', j-1, jb,
534 $ n-k, -one, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
535 $ one, a( 1, j ), lda )
536 50 CONTINUE
537*
538* Put U12 in standard form by partially undoing the interchanges
539* in columns k+1:n
540*
541 j = k + 1
542 60 CONTINUE
543*
544 kstep = 1
545 jp1 = 1
546 jj = j
547 jp2 = ipiv( j )
548 IF( jp2.LT.0 ) THEN
549 jp2 = -jp2
550 j = j + 1
551 jp1 = -ipiv( j )
552 kstep = 2
553 END IF
554*
555 j = j + 1
556 IF( jp2.NE.jj .AND. j.LE.n )
557 $ CALL sswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
558 jj = j - 1
559 IF( jp1.NE.jj .AND. kstep.EQ.2 )
560 $ CALL sswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
561 IF( j.LE.n )
562 $ GO TO 60
563*
564* Set KB to the number of columns factorized
565*
566 kb = n - k
567*
568 ELSE
569*
570* Factorize the leading columns of A using the lower triangle
571* of A and working forwards, and compute the matrix W = L21*D
572* for use in updating A22
573*
574* K is the main loop index, increasing from 1 in steps of 1 or 2
575*
576 k = 1
577 70 CONTINUE
578*
579* Exit from loop
580*
581 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
582 $ GO TO 90
583*
584 kstep = 1
585 p = k
586*
587* Copy column K of A to column K of W and update it
588*
589 CALL scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
590 IF( k.GT.1 )
591 $ CALL sgemv( 'No transpose', n-k+1, k-1, -one, a( k, 1 ),
592 $ lda, w( k, 1 ), ldw, one, w( k, k ), 1 )
593*
594* Determine rows and columns to be interchanged and whether
595* a 1-by-1 or 2-by-2 pivot block will be used
596*
597 absakk = abs( w( k, k ) )
598*
599* IMAX is the row-index of the largest off-diagonal element in
600* column K, and COLMAX is its absolute value.
601* Determine both COLMAX and IMAX.
602*
603 IF( k.LT.n ) THEN
604 imax = k + isamax( n-k, w( k+1, k ), 1 )
605 colmax = abs( w( imax, k ) )
606 ELSE
607 colmax = zero
608 END IF
609*
610 IF( max( absakk, colmax ).EQ.zero ) THEN
611*
612* Column K is zero or underflow: set INFO and continue
613*
614 IF( info.EQ.0 )
615 $ info = k
616 kp = k
617 CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
618 ELSE
619*
620* ============================================================
621*
622* Test for interchange
623*
624* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
625* (used to handle NaN and Inf)
626*
627 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
628*
629* no interchange, use 1-by-1 pivot block
630*
631 kp = k
632*
633 ELSE
634*
635 done = .false.
636*
637* Loop until pivot found
638*
639 72 CONTINUE
640*
641* Begin pivot search loop body
642*
643*
644* Copy column IMAX to column K+1 of W and update it
645*
646 CALL scopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
647 CALL scopy( n-imax+1, a( imax, imax ), 1,
648 $ w( imax, k+1 ), 1 )
649 IF( k.GT.1 )
650 $ CALL sgemv( 'No transpose', n-k+1, k-1, -one,
651 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
652 $ one, w( k, k+1 ), 1 )
653*
654* JMAX is the column-index of the largest off-diagonal
655* element in row IMAX, and ROWMAX is its absolute value.
656* Determine both ROWMAX and JMAX.
657*
658 IF( imax.NE.k ) THEN
659 jmax = k - 1 + isamax( imax-k, w( k, k+1 ), 1 )
660 rowmax = abs( w( jmax, k+1 ) )
661 ELSE
662 rowmax = zero
663 END IF
664*
665 IF( imax.LT.n ) THEN
666 itemp = imax + isamax( n-imax, w( imax+1, k+1 ), 1)
667 stemp = abs( w( itemp, k+1 ) )
668 IF( stemp.GT.rowmax ) THEN
669 rowmax = stemp
670 jmax = itemp
671 END IF
672 END IF
673*
674* Equivalent to testing for
675* ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
676* (used to handle NaN and Inf)
677*
678 IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
679 $ THEN
680*
681* interchange rows and columns K and IMAX,
682* use 1-by-1 pivot block
683*
684 kp = imax
685*
686* copy column K+1 of W to column K of W
687*
688 CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
689*
690 done = .true.
691*
692* Equivalent to testing for ROWMAX.EQ.COLMAX,
693* (used to handle NaN and Inf)
694*
695 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
696 $ THEN
697*
698* interchange rows and columns K+1 and IMAX,
699* use 2-by-2 pivot block
700*
701 kp = imax
702 kstep = 2
703 done = .true.
704 ELSE
705*
706* Pivot not found: set params and repeat
707*
708 p = imax
709 colmax = rowmax
710 imax = jmax
711*
712* Copy updated JMAXth (next IMAXth) column to Kth of W
713*
714 CALL scopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
715*
716 END IF
717*
718* End pivot search loop body
719*
720 IF( .NOT. done ) GOTO 72
721*
722 END IF
723*
724* ============================================================
725*
726 kk = k + kstep - 1
727*
728 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
729*
730* Copy non-updated column K to column P
731*
732 CALL scopy( p-k, a( k, k ), 1, a( p, k ), lda )
733 CALL scopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
734*
735* Interchange rows K and P in first K columns of A
736* and first K+1 columns of W
737*
738 CALL sswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
739 CALL sswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
740 END IF
741*
742* Updated column KP is already stored in column KK of W
743*
744 IF( kp.NE.kk ) THEN
745*
746* Copy non-updated column KK to column KP
747*
748 a( kp, k ) = a( kk, k )
749 CALL scopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
750 CALL scopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
751*
752* Interchange rows KK and KP in first KK columns of A and W
753*
754 CALL sswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
755 CALL sswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
756 END IF
757*
758 IF( kstep.EQ.1 ) THEN
759*
760* 1-by-1 pivot block D(k): column k of W now holds
761*
762* W(k) = L(k)*D(k)
763*
764* where L(k) is the k-th column of L
765*
766* Store L(k) in column k of A
767*
768 CALL scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
769 IF( k.LT.n ) THEN
770 IF( abs( a( k, k ) ).GE.sfmin ) THEN
771 r1 = one / a( k, k )
772 CALL sscal( n-k, r1, a( k+1, k ), 1 )
773 ELSE IF( a( k, k ).NE.zero ) THEN
774 DO 74 ii = k + 1, n
775 a( ii, k ) = a( ii, k ) / a( k, k )
776 74 CONTINUE
777 END IF
778 END IF
779*
780 ELSE
781*
782* 2-by-2 pivot block D(k): columns k and k+1 of W now hold
783*
784* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
785*
786* where L(k) and L(k+1) are the k-th and (k+1)-th columns
787* of L
788*
789 IF( k.LT.n-1 ) THEN
790*
791* Store L(k) and L(k+1) in columns k and k+1 of A
792*
793 d21 = w( k+1, k )
794 d11 = w( k+1, k+1 ) / d21
795 d22 = w( k, k ) / d21
796 t = one / ( d11*d22-one )
797 DO 80 j = k + 2, n
798 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
799 $ d21 )
800 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
801 $ d21 )
802 80 CONTINUE
803 END IF
804*
805* Copy D(k) to A
806*
807 a( k, k ) = w( k, k )
808 a( k+1, k ) = w( k+1, k )
809 a( k+1, k+1 ) = w( k+1, k+1 )
810 END IF
811 END IF
812*
813* Store details of the interchanges in IPIV
814*
815 IF( kstep.EQ.1 ) THEN
816 ipiv( k ) = kp
817 ELSE
818 ipiv( k ) = -p
819 ipiv( k+1 ) = -kp
820 END IF
821*
822* Increase K and return to the start of the main loop
823*
824 k = k + kstep
825 GO TO 70
826*
827 90 CONTINUE
828*
829* Update the lower triangle of A22 (= A(k:n,k:n)) as
830*
831* A22 := A22 - L21*D*L21**T = A22 - L21*W**T
832*
833* computing blocks of NB columns at a time
834*
835 DO 110 j = k, n, nb
836 jb = min( nb, n-j+1 )
837*
838* Update the lower triangle of the diagonal block
839*
840 DO 100 jj = j, j + jb - 1
841 CALL sgemv( 'No transpose', j+jb-jj, k-1, -one,
842 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
843 $ a( jj, jj ), 1 )
844 100 CONTINUE
845*
846* Update the rectangular subdiagonal block
847*
848 IF( j+jb.LE.n )
849 $ CALL sgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
850 $ k-1, -one, a( j+jb, 1 ), lda, w( j, 1 ), ldw,
851 $ one, a( j+jb, j ), lda )
852 110 CONTINUE
853*
854* Put L21 in standard form by partially undoing the interchanges
855* in columns 1:k-1
856*
857 j = k - 1
858 120 CONTINUE
859*
860 kstep = 1
861 jp1 = 1
862 jj = j
863 jp2 = ipiv( j )
864 IF( jp2.LT.0 ) THEN
865 jp2 = -jp2
866 j = j - 1
867 jp1 = -ipiv( j )
868 kstep = 2
869 END IF
870*
871 j = j - 1
872 IF( jp2.NE.jj .AND. j.GE.1 )
873 $ CALL sswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
874 jj = j + 1
875 IF( jp1.NE.jj .AND. kstep.EQ.2 )
876 $ CALL sswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
877 IF( j.GE.1 )
878 $ GO TO 120
879*
880* Set KB to the number of columns factorized
881*
882 kb = k - 1
883*
884 END IF
885 RETURN
886*
887* End of SLASYF_ROOK
888*
889 END
subroutine slasyf_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
SLASYF_ROOK computes a partial factorization of a real symmetric matrix using the bounded Bunch-Kaufm...
Definition: slasyf_rook.f:184
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 sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187