LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsytf2_rook.f
Go to the documentation of this file.
1*> \brief \b DSYTF2_ROOK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method (unblocked algorithm).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSYTF2_ROOK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytf2_rook.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytf2_rook.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytf2_rook.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, LDA, N
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* DOUBLE PRECISION A( LDA, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DSYTF2_ROOK computes the factorization of a real symmetric matrix A
39*> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method:
40*>
41*> A = U*D*U**T or A = L*D*L**T
42*>
43*> where U (or L) is a product of permutation and unit upper (lower)
44*> triangular matrices, U**T is the transpose of U, and D is symmetric and
45*> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
46*>
47*> This is the unblocked version of the algorithm, calling Level 2 BLAS.
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] UPLO
54*> \verbatim
55*> UPLO is CHARACTER*1
56*> Specifies whether the upper or lower triangular part of the
57*> symmetric matrix A is stored:
58*> = 'U': Upper triangular
59*> = 'L': Lower triangular
60*> \endverbatim
61*>
62*> \param[in] N
63*> \verbatim
64*> N is INTEGER
65*> The order of the matrix A. N >= 0.
66*> \endverbatim
67*>
68*> \param[in,out] A
69*> \verbatim
70*> A is DOUBLE PRECISION array, dimension (LDA,N)
71*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
72*> n-by-n upper triangular part of A contains the upper
73*> triangular part of the matrix A, and the strictly lower
74*> triangular part of A is not referenced. If UPLO = 'L', the
75*> leading n-by-n lower triangular part of A contains the lower
76*> triangular part of the matrix A, and the strictly upper
77*> triangular part of A is not referenced.
78*>
79*> On exit, the block diagonal matrix D and the multipliers used
80*> to obtain the factor U or L (see below for further details).
81*> \endverbatim
82*>
83*> \param[in] LDA
84*> \verbatim
85*> LDA is INTEGER
86*> The leading dimension of the array A. LDA >= max(1,N).
87*> \endverbatim
88*>
89*> \param[out] IPIV
90*> \verbatim
91*> IPIV is INTEGER array, dimension (N)
92*> Details of the interchanges and the block structure of D.
93*>
94*> If UPLO = 'U':
95*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
96*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
97*>
98*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
99*> columns k and -IPIV(k) were interchanged and rows and
100*> columns k-1 and -IPIV(k-1) were inerchaged,
101*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
102*>
103*> If UPLO = 'L':
104*> If IPIV(k) > 0, then rows and columns k and IPIV(k)
105*> were interchanged and D(k,k) is a 1-by-1 diagonal block.
106*>
107*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
108*> columns k and -IPIV(k) were interchanged and rows and
109*> columns k+1 and -IPIV(k+1) were inerchaged,
110*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
111*> \endverbatim
112*>
113*> \param[out] INFO
114*> \verbatim
115*> INFO is INTEGER
116*> = 0: successful exit
117*> < 0: if INFO = -k, the k-th argument had an illegal value
118*> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
119*> has been completed, but the block diagonal matrix D is
120*> exactly singular, and division by zero will occur if it
121*> is used to solve a system of equations.
122*> \endverbatim
123*
124* Authors:
125* ========
126*
127*> \author Univ. of Tennessee
128*> \author Univ. of California Berkeley
129*> \author Univ. of Colorado Denver
130*> \author NAG Ltd.
131*
132*> \ingroup hetf2_rook
133*
134*> \par Further Details:
135* =====================
136*>
137*> \verbatim
138*>
139*> If UPLO = 'U', then A = U*D*U**T, where
140*> U = P(n)*U(n)* ... *P(k)U(k)* ...,
141*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
142*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
143*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
144*> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
145*> that if the diagonal block D(k) is of order s (s = 1 or 2), then
146*>
147*> ( I v 0 ) k-s
148*> U(k) = ( 0 I 0 ) s
149*> ( 0 0 I ) n-k
150*> k-s s n-k
151*>
152*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
153*> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
154*> and A(k,k), and v overwrites A(1:k-2,k-1:k).
155*>
156*> If UPLO = 'L', then A = L*D*L**T, where
157*> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
158*> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
159*> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
160*> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
161*> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
162*> that if the diagonal block D(k) is of order s (s = 1 or 2), then
163*>
164*> ( I 0 0 ) k-1
165*> L(k) = ( 0 I 0 ) s
166*> ( 0 v I ) n-k-s+1
167*> k-1 s n-k-s+1
168*>
169*> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
170*> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
171*> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
172*> \endverbatim
173*
174*> \par Contributors:
175* ==================
176*>
177*> \verbatim
178*>
179*> November 2013, Igor Kozachenko,
180*> Computer Science Division,
181*> University of California, Berkeley
182*>
183*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
184*> School of Mathematics,
185*> University of Manchester
186*>
187*> 01-01-96 - Based on modifications by
188*> J. Lewis, Boeing Computer Services Company
189*> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville abd , USA
190*> \endverbatim
191*
192* =====================================================================
193 SUBROUTINE dsytf2_rook( UPLO, N, A, LDA, IPIV, INFO )
194*
195* -- LAPACK computational routine --
196* -- LAPACK is a software package provided by Univ. of Tennessee, --
197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198*
199* .. Scalar Arguments ..
200 CHARACTER UPLO
201 INTEGER INFO, LDA, N
202* ..
203* .. Array Arguments ..
204 INTEGER IPIV( * )
205 DOUBLE PRECISION A( LDA, * )
206* ..
207*
208* =====================================================================
209*
210* .. Parameters ..
211 DOUBLE PRECISION ZERO, ONE
212 parameter( zero = 0.0d+0, one = 1.0d+0 )
213 DOUBLE PRECISION EIGHT, SEVTEN
214 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
215* ..
216* .. Local Scalars ..
217 LOGICAL UPPER, DONE
218 INTEGER I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
219 $ P, II
220 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22,
221 $ ROWMAX, DTEMP, T, WK, WKM1, WKP1, SFMIN
222* ..
223* .. External Functions ..
224 LOGICAL LSAME
225 INTEGER IDAMAX
226 DOUBLE PRECISION DLAMCH
227 EXTERNAL lsame, idamax, dlamch
228* ..
229* .. External Subroutines ..
230 EXTERNAL dscal, dswap, dsyr, xerbla
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC abs, max, sqrt
234* ..
235* .. Executable Statements ..
236*
237* Test the input parameters.
238*
239 info = 0
240 upper = lsame( uplo, 'U' )
241 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
242 info = -1
243 ELSE IF( n.LT.0 ) THEN
244 info = -2
245 ELSE IF( lda.LT.max( 1, n ) ) THEN
246 info = -4
247 END IF
248 IF( info.NE.0 ) THEN
249 CALL xerbla( 'DSYTF2_ROOK', -info )
250 RETURN
251 END IF
252*
253* Initialize ALPHA for use in choosing pivot block size.
254*
255 alpha = ( one+sqrt( sevten ) ) / eight
256*
257* Compute machine safe minimum
258*
259 sfmin = dlamch( 'S' )
260*
261 IF( upper ) THEN
262*
263* Factorize A as U*D*U**T using the upper triangle of A
264*
265* K is the main loop index, decreasing from N to 1 in steps of
266* 1 or 2
267*
268 k = n
269 10 CONTINUE
270*
271* If K < 1, exit from loop
272*
273 IF( k.LT.1 )
274 $ GO TO 70
275 kstep = 1
276 p = k
277*
278* Determine rows and columns to be interchanged and whether
279* a 1-by-1 or 2-by-2 pivot block will be used
280*
281 absakk = abs( a( k, k ) )
282*
283* IMAX is the row-index of the largest off-diagonal element in
284* column K, and COLMAX is its absolute value.
285* Determine both COLMAX and IMAX.
286*
287 IF( k.GT.1 ) THEN
288 imax = idamax( k-1, a( 1, k ), 1 )
289 colmax = abs( a( imax, k ) )
290 ELSE
291 colmax = zero
292 END IF
293*
294 IF( (max( absakk, colmax ).EQ.zero) ) THEN
295*
296* Column K is zero or underflow: set INFO and continue
297*
298 IF( info.EQ.0 )
299 $ info = k
300 kp = k
301 ELSE
302*
303* Test for interchange
304*
305* Equivalent to testing for (used to handle NaN and Inf)
306* ABSAKK.GE.ALPHA*COLMAX
307*
308 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
309*
310* no interchange,
311* use 1-by-1 pivot block
312*
313 kp = k
314 ELSE
315*
316 done = .false.
317*
318* Loop until pivot found
319*
320 12 CONTINUE
321*
322* Begin pivot search loop body
323*
324* JMAX is the column-index of the largest off-diagonal
325* element in row IMAX, and ROWMAX is its absolute value.
326* Determine both ROWMAX and JMAX.
327*
328 IF( imax.NE.k ) THEN
329 jmax = imax + idamax( k-imax, a( imax, imax+1 ),
330 $ lda )
331 rowmax = abs( a( imax, jmax ) )
332 ELSE
333 rowmax = zero
334 END IF
335*
336 IF( imax.GT.1 ) THEN
337 itemp = idamax( imax-1, a( 1, imax ), 1 )
338 dtemp = abs( a( itemp, imax ) )
339 IF( dtemp.GT.rowmax ) THEN
340 rowmax = dtemp
341 jmax = itemp
342 END IF
343 END IF
344*
345* Equivalent to testing for (used to handle NaN and Inf)
346* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
347*
348 IF( .NOT.( abs( a( imax, imax ) ).LT.alpha*rowmax ) )
349 $ THEN
350*
351* interchange rows and columns K and IMAX,
352* use 1-by-1 pivot block
353*
354 kp = imax
355 done = .true.
356*
357* Equivalent to testing for ROWMAX .EQ. COLMAX,
358* used to handle NaN and Inf
359*
360 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
361*
362* interchange rows and columns K+1 and IMAX,
363* use 2-by-2 pivot block
364*
365 kp = imax
366 kstep = 2
367 done = .true.
368 ELSE
369*
370* Pivot NOT found, set variables and repeat
371*
372 p = imax
373 colmax = rowmax
374 imax = jmax
375 END IF
376*
377* End pivot search loop body
378*
379 IF( .NOT. done ) GOTO 12
380*
381 END IF
382*
383* Swap TWO rows and TWO columns
384*
385* First swap
386*
387 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
388*
389* Interchange rows and column K and P in the leading
390* submatrix A(1:k,1:k) if we have a 2-by-2 pivot
391*
392 IF( p.GT.1 )
393 $ CALL dswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
394 IF( p.LT.(k-1) )
395 $ CALL dswap( k-p-1, a( p+1, k ), 1, a( p, p+1 ),
396 $ lda )
397 t = a( k, k )
398 a( k, k ) = a( p, p )
399 a( p, p ) = t
400 END IF
401*
402* Second swap
403*
404 kk = k - kstep + 1
405 IF( kp.NE.kk ) THEN
406*
407* Interchange rows and columns KK and KP in the leading
408* submatrix A(1:k,1:k)
409*
410 IF( kp.GT.1 )
411 $ CALL dswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
412 IF( ( kk.GT.1 ) .AND. ( kp.LT.(kk-1) ) )
413 $ CALL dswap( kk-kp-1, a( kp+1, kk ), 1, a( kp, kp+1 ),
414 $ lda )
415 t = a( kk, kk )
416 a( kk, kk ) = a( kp, kp )
417 a( kp, kp ) = t
418 IF( kstep.EQ.2 ) THEN
419 t = a( k-1, k )
420 a( k-1, k ) = a( kp, k )
421 a( kp, k ) = t
422 END IF
423 END IF
424*
425* Update the leading submatrix
426*
427 IF( kstep.EQ.1 ) THEN
428*
429* 1-by-1 pivot block D(k): column k now holds
430*
431* W(k) = U(k)*D(k)
432*
433* where U(k) is the k-th column of U
434*
435 IF( k.GT.1 ) THEN
436*
437* Perform a rank-1 update of A(1:k-1,1:k-1) and
438* store U(k) in column k
439*
440 IF( abs( a( k, k ) ).GE.sfmin ) THEN
441*
442* Perform a rank-1 update of A(1:k-1,1:k-1) as
443* A := A - U(k)*D(k)*U(k)**T
444* = A - W(k)*1/D(k)*W(k)**T
445*
446 d11 = one / a( k, k )
447 CALL dsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
448*
449* Store U(k) in column k
450*
451 CALL dscal( k-1, d11, a( 1, k ), 1 )
452 ELSE
453*
454* Store L(k) in column K
455*
456 d11 = a( k, k )
457 DO 16 ii = 1, k - 1
458 a( ii, k ) = a( ii, k ) / d11
459 16 CONTINUE
460*
461* Perform a rank-1 update of A(k+1:n,k+1:n) as
462* A := A - U(k)*D(k)*U(k)**T
463* = A - W(k)*(1/D(k))*W(k)**T
464* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
465*
466 CALL dsyr( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
467 END IF
468 END IF
469*
470 ELSE
471*
472* 2-by-2 pivot block D(k): columns k and k-1 now hold
473*
474* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
475*
476* where U(k) and U(k-1) are the k-th and (k-1)-th columns
477* of U
478*
479* Perform a rank-2 update of A(1:k-2,1:k-2) as
480*
481* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
482* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
483*
484* and store L(k) and L(k+1) in columns k and k+1
485*
486 IF( k.GT.2 ) THEN
487*
488 d12 = a( k-1, k )
489 d22 = a( k-1, k-1 ) / d12
490 d11 = a( k, k ) / d12
491 t = one / ( d11*d22-one )
492*
493 DO 30 j = k - 2, 1, -1
494*
495 wkm1 = t*( d11*a( j, k-1 )-a( j, k ) )
496 wk = t*( d22*a( j, k )-a( j, k-1 ) )
497*
498 DO 20 i = j, 1, -1
499 a( i, j ) = a( i, j ) - (a( i, k ) / d12 )*wk -
500 $ ( a( i, k-1 ) / d12 )*wkm1
501 20 CONTINUE
502*
503* Store U(k) and U(k-1) in cols k and k-1 for row J
504*
505 a( j, k ) = wk / d12
506 a( j, k-1 ) = wkm1 / d12
507*
508 30 CONTINUE
509*
510 END IF
511*
512 END IF
513 END IF
514*
515* Store details of the interchanges in IPIV
516*
517 IF( kstep.EQ.1 ) THEN
518 ipiv( k ) = kp
519 ELSE
520 ipiv( k ) = -p
521 ipiv( k-1 ) = -kp
522 END IF
523*
524* Decrease K and return to the start of the main loop
525*
526 k = k - kstep
527 GO TO 10
528*
529 ELSE
530*
531* Factorize A as L*D*L**T using the lower triangle of A
532*
533* K is the main loop index, increasing from 1 to N in steps of
534* 1 or 2
535*
536 k = 1
537 40 CONTINUE
538*
539* If K > N, exit from loop
540*
541 IF( k.GT.n )
542 $ GO TO 70
543 kstep = 1
544 p = k
545*
546* Determine rows and columns to be interchanged and whether
547* a 1-by-1 or 2-by-2 pivot block will be used
548*
549 absakk = abs( a( k, k ) )
550*
551* IMAX is the row-index of the largest off-diagonal element in
552* column K, and COLMAX is its absolute value.
553* Determine both COLMAX and IMAX.
554*
555 IF( k.LT.n ) THEN
556 imax = k + idamax( n-k, a( k+1, k ), 1 )
557 colmax = abs( a( imax, k ) )
558 ELSE
559 colmax = zero
560 END IF
561*
562 IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
563*
564* Column K is zero or underflow: set INFO and continue
565*
566 IF( info.EQ.0 )
567 $ info = k
568 kp = k
569 ELSE
570*
571* Test for interchange
572*
573* Equivalent to testing for (used to handle NaN and Inf)
574* ABSAKK.GE.ALPHA*COLMAX
575*
576 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
577*
578* no interchange, use 1-by-1 pivot block
579*
580 kp = k
581 ELSE
582*
583 done = .false.
584*
585* Loop until pivot found
586*
587 42 CONTINUE
588*
589* Begin pivot search loop body
590*
591* JMAX is the column-index of the largest off-diagonal
592* element in row IMAX, and ROWMAX is its absolute value.
593* Determine both ROWMAX and JMAX.
594*
595 IF( imax.NE.k ) THEN
596 jmax = k - 1 + idamax( imax-k, a( imax, k ), lda )
597 rowmax = abs( a( imax, jmax ) )
598 ELSE
599 rowmax = zero
600 END IF
601*
602 IF( imax.LT.n ) THEN
603 itemp = imax + idamax( n-imax, a( imax+1, imax ),
604 $ 1 )
605 dtemp = abs( a( itemp, imax ) )
606 IF( dtemp.GT.rowmax ) THEN
607 rowmax = dtemp
608 jmax = itemp
609 END IF
610 END IF
611*
612* Equivalent to testing for (used to handle NaN and Inf)
613* ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
614*
615 IF( .NOT.( abs( a( imax, imax ) ).LT.alpha*rowmax ) )
616 $ THEN
617*
618* interchange rows and columns K and IMAX,
619* use 1-by-1 pivot block
620*
621 kp = imax
622 done = .true.
623*
624* Equivalent to testing for ROWMAX .EQ. COLMAX,
625* used to handle NaN and Inf
626*
627 ELSE IF( ( p.EQ.jmax ).OR.( rowmax.LE.colmax ) ) THEN
628*
629* interchange rows and columns K+1 and IMAX,
630* use 2-by-2 pivot block
631*
632 kp = imax
633 kstep = 2
634 done = .true.
635 ELSE
636*
637* Pivot NOT found, set variables and repeat
638*
639 p = imax
640 colmax = rowmax
641 imax = jmax
642 END IF
643*
644* End pivot search loop body
645*
646 IF( .NOT. done ) GOTO 42
647*
648 END IF
649*
650* Swap TWO rows and TWO columns
651*
652* First swap
653*
654 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
655*
656* Interchange rows and column K and P in the trailing
657* submatrix A(k:n,k:n) if we have a 2-by-2 pivot
658*
659 IF( p.LT.n )
660 $ CALL dswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
661 IF( p.GT.(k+1) )
662 $ CALL dswap( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
663 t = a( k, k )
664 a( k, k ) = a( p, p )
665 a( p, p ) = t
666 END IF
667*
668* Second swap
669*
670 kk = k + kstep - 1
671 IF( kp.NE.kk ) THEN
672*
673* Interchange rows and columns KK and KP in the trailing
674* submatrix A(k:n,k:n)
675*
676 IF( kp.LT.n )
677 $ CALL dswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
678 IF( ( kk.LT.n ) .AND. ( kp.GT.(kk+1) ) )
679 $ CALL dswap( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
680 $ lda )
681 t = a( kk, kk )
682 a( kk, kk ) = a( kp, kp )
683 a( kp, kp ) = t
684 IF( kstep.EQ.2 ) THEN
685 t = a( k+1, k )
686 a( k+1, k ) = a( kp, k )
687 a( kp, k ) = t
688 END IF
689 END IF
690*
691* Update the trailing submatrix
692*
693 IF( kstep.EQ.1 ) THEN
694*
695* 1-by-1 pivot block D(k): column k now holds
696*
697* W(k) = L(k)*D(k)
698*
699* where L(k) is the k-th column of L
700*
701 IF( k.LT.n ) THEN
702*
703* Perform a rank-1 update of A(k+1:n,k+1:n) and
704* store L(k) in column k
705*
706 IF( abs( a( k, k ) ).GE.sfmin ) THEN
707*
708* Perform a rank-1 update of A(k+1:n,k+1:n) as
709* A := A - L(k)*D(k)*L(k)**T
710* = A - W(k)*(1/D(k))*W(k)**T
711*
712 d11 = one / a( k, k )
713 CALL dsyr( uplo, n-k, -d11, a( k+1, k ), 1,
714 $ a( k+1, k+1 ), lda )
715*
716* Store L(k) in column k
717*
718 CALL dscal( n-k, d11, a( k+1, k ), 1 )
719 ELSE
720*
721* Store L(k) in column k
722*
723 d11 = a( k, k )
724 DO 46 ii = k + 1, n
725 a( ii, k ) = a( ii, k ) / d11
726 46 CONTINUE
727*
728* Perform a rank-1 update of A(k+1:n,k+1:n) as
729* A := A - L(k)*D(k)*L(k)**T
730* = A - W(k)*(1/D(k))*W(k)**T
731* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
732*
733 CALL dsyr( uplo, n-k, -d11, a( k+1, k ), 1,
734 $ a( k+1, k+1 ), lda )
735 END IF
736 END IF
737*
738 ELSE
739*
740* 2-by-2 pivot block D(k): columns k and k+1 now hold
741*
742* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
743*
744* where L(k) and L(k+1) are the k-th and (k+1)-th columns
745* of L
746*
747*
748* Perform a rank-2 update of A(k+2:n,k+2:n) as
749*
750* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
751* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
752*
753* and store L(k) and L(k+1) in columns k and k+1
754*
755 IF( k.LT.n-1 ) THEN
756*
757 d21 = a( k+1, k )
758 d11 = a( k+1, k+1 ) / d21
759 d22 = a( k, k ) / d21
760 t = one / ( d11*d22-one )
761*
762 DO 60 j = k + 2, n
763*
764* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
765*
766 wk = t*( d11*a( j, k )-a( j, k+1 ) )
767 wkp1 = t*( d22*a( j, k+1 )-a( j, k ) )
768*
769* Perform a rank-2 update of A(k+2:n,k+2:n)
770*
771 DO 50 i = j, n
772 a( i, j ) = a( i, j ) - ( a( i, k ) / d21 )*wk -
773 $ ( a( i, k+1 ) / d21 )*wkp1
774 50 CONTINUE
775*
776* Store L(k) and L(k+1) in cols k and k+1 for row J
777*
778 a( j, k ) = wk / d21
779 a( j, k+1 ) = wkp1 / d21
780*
781 60 CONTINUE
782*
783 END IF
784*
785 END IF
786 END IF
787*
788* Store details of the interchanges in IPIV
789*
790 IF( kstep.EQ.1 ) THEN
791 ipiv( k ) = kp
792 ELSE
793 ipiv( k ) = -p
794 ipiv( k+1 ) = -kp
795 END IF
796*
797* Increase K and return to the start of the main loop
798*
799 k = k + kstep
800 GO TO 40
801*
802 END IF
803*
804 70 CONTINUE
805*
806 RETURN
807*
808* End of DSYTF2_ROOK
809*
810 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsyr(uplo, n, alpha, x, incx, a, lda)
DSYR
Definition dsyr.f:132
subroutine dsytf2_rook(uplo, n, a, lda, ipiv, info)
DSYTF2_ROOK computes the factorization of a real symmetric indefinite matrix using the bounded Bunch-...
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82